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Abstract 


The  purpose  of  this  research  was  to  develop  a  highly  accurate 
computational  method  for  calculating  the  nonequilibrium  radiative  heat  transfer 
within  reentry  shock  layers.  The  nonequilibrium  state  of  the  flowfield  was 
obtained  by  using  the  multispecies  multitemperature  nonequilibrium  flow  solver 
NH7AIR  which  is  capable  of  separately  tracking  the  vibrational  energy  of  each 
diatomic  species  and  the  energy  of  the  free  and  bound  electrons.  The  calculation 
of  radiative  heat  transfer  was  performed  by  utilizing  the  detailed  line-by-line 
spectral  radiation  solver  SPRADIAN.  Two  radiative  transport  schemes  were 
implemented  in  this  coupled  code.  The  first  scheme  was  based  on  a 
straightforward  application  of  the  standard  tangent  slab  solution  method.  The 
second  scheme  was  based  on  the  conservation  of  radiative  energy  and  resulted  in 
a  finite  volume  method  for  radiative  heat  transfer  (FVMR).  Data  from  the  FIRE 
11  flight  experiment  were  used  to  validate  the  coupled  radiation-gasdynamic 
solver.  Coupled  results  exhibited  a  high  degree  of  agreement  with  experimental 
data.  The  utility  of  the  FVMR  scheme  was  also  examined  in  an  uncoupled 
implementation  and  shows  promise  for  future  implementation  in  a  coupled 
setting.  Together,  the  enhancement  of  the  nonequilibrium  thermal  modeling,  the 
use  of  a  highly  accurate  spectral  radiation  solver  and  the  development  of  a 
conservative  scheme  for  radiative  transport  constitute  a  significant  improvement 
in  current  capabilities  available  for  modeling  the  radiating  nonequilibrium  shock 
layers  which  accompany  reentry  flight. 
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Preface 


Perseverance  is  more  prevailing  than  violence;  and  many  things  which  cannot  be  overcome  when  they  are 
together,  yield  themselves  up  when  taken  little  by  little. 

-Plutarch,  Parallel  Lives,  Sertorius,  1st  c. 

Analysis — definition:  noun.  1 .  a  detailed  examination  of  the  elements  or  structure  of  something,  2.  the 
separation  of  something  into  its  constituent  elements 

-Oxford  English  Dictionary,  3rd  ed. ,  revised,  2008 

Plutarch‘s  observation  regarding  the  career  of  the  Roman  general  Quintus 

Sertorius  holds  fairly  well  by  analogy  to  the  task  of  analysis  set  before  the  student — and 
perhaps  even  especially  so  for  anyone  conducting  research  within  the  field  of 
hypersonics.  The  implications  of  hypersonic  flight  for  mankind  are  indeed  exciting,  and 
many  large,  national  programs  have  been  enthusiastically  advanced  in  pursuit  of  this  lofty 
goal.  However,  many  of  them  have  failed  to  persevere  under  various  technological, 
institutional,  fiscal  and  political  constraints;  that  this  is  evident  may  be  seen  from  the  fits- 
and-starts  history  of  hypersonics  to  date.  The  great  challenge  is  that  hypersonic  flows  are 
intrinsically  multiphysical  in  nature  and  encompasses  numerous  physical  phenomena 
issuing  from  many  distinct  fields  of  inquiry:  gasdynamics,  gaskinetic  theory, 
thermodynamics,  turbulence,  material  science,  chemistry,  radiative  transport,  quantum 
mechanics,  et  cetera.  Thus,  those  who  would  master  the  theoretical  and  technological 
challenges  of  the  field  must,  like  the  Roman  generals  who  shrewdly  expanded  upon  the 
earth  a  vast  and  disparate  empire,  carefully  examine  in  detail  these  varied  phenomena 
with  an  eye  to  their  characteristics  individually,  as  well  as  their  characteristics  in 
interaction  with  one  another.  Such  is  the  nature  of  even  modest  applications  of 
theoretical  knowledge  to  the  design  and  analysis  of  real-world  hypersonic  flight 
technologies. 
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COUPLED  RADIATION -GASDYNAMIC  SOLUTION  METHOD  FOR 
HYPERSONIC  SHOCK  LAYERS  IN  THERMOCHEMICAL  NONEQUILIBRIUM 

I.  Introduction 

At  the  time  of  this  writing,  the  fact  that  national  interest  in  hypersonics  is  alive 
and  well — or  at  least  making  one  of  its  periodic  comebacks — is  evident  by  the  hundred- 
million-dollar  expenditures  which  the  DoD  has  made  in  this  past  decade  to  reinvigorate 
our  national  competency  in  hypersonics  (OSD,  2008).  Recent  comments  from  Dr.  Dahm, 
former  Chief  Scientist  of  the  Air  Force,  are  a  further  indication  of  this  upward  swing  in 
interest.  Dr.  Dahm  has  been  a  potent  advocate  for  hypersonics  within  the  DoD  and  has 
labored  to  aid  decision  makers  in  understanding  the  contemporary  implications  of 
hypersonic  systems,  noting  that  —  .this  wouldn‘t  just  do  what  we  do  today  faster.  We 
could  do  things  differently”  (Barnes,  2010).  In  other  words,  the  fielding  of  operational 
hypersonics  systems  by  the  United  States  or  a  competitor  nation  would  constitute  a 
disruptive  paradigm  shift  (Borger,  2007),  substantially  affecting  the  way  in  which  wars 
are  prosecuted  and  the  homeland  defended.  However,  Dr.  Dahm  and  other  proponents  of 
the  military  utility  of  such  hypersonic  systems  must  wait  patiently,  if  not  a  bit  anxiously, 
for  the  completion  of  the  requisite  basic  research  and  technology  development.  Perhaps 
this  anxiety  is  understandable,  since  the  corporate  memory  of  the  hypersonics  community 
is  haunted  by  the  ghosts  of  many  canceled  programs  which  either  failed  to  perform 
according  to  stakeholder  expectations  or  which  were  otherwise  deemed  too  risky  to 
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pursue  (Heppenheimer,  2007;  Hallion  2005;  Anderson,  1984).  Yet,  despite  towering 
technical  challenges  and  many  perceived  setbacks,  the  Air  Force  has  throughout  its 
history  maintained  some  vision  for  the  utilization  of  hypersonic  flight  vehicles,  although 
reformulated  in  various  ways  according  to  the  perceived  needs  of  the  day.  In  recent 
years,  these  needs  have  come  to  be  redefined  within  the  purview  of  the  Operationally 
Responsive  Space  and  Global  Strike/Global  Persistent  Attack  CONOPS,  for  which  some 
sort  of  reusable,  air-breathing  or  rocket-launched  platform  is  typically  envisioned  (Fuchs, 
et  al.,  2000;  Tichkoff,  et  al.,  1998;  McCall,  et  al.,  1995;  McLucas,  et  ah,  1989). 

Whatever  vehicle  concept  is  chosen  to  enable  these  CONOPS,  the  basic  design 
and  analysis  tasks  in  the  development  of  a  hypersonic  system  remain  fundamentally 
unchanged.  This  fact  is  quite  evident  from  even  a  cursory  reading  of  Heppenheimer‘s 
Facing  the  Heat  Barrier:  a  History  of  Hypersonics.  Whether  discussing  X-15,  Dyna-Soar, 
Apollo,  NASP,  or  X-51,  the  same  technical  problem  areas  are  addressed  again  and  again, 
namely:  propulsion,  materials,  structures,  transition/turbulence,  control,  and,  finally, 
thermal  management. 

The  scope  of  the  accomplished  research  is  first  within  the  bounds  of  the  latter 
problem  area.  As  Heppenheimer4  s  selected  title  indicates,  one  of  the  most  significant 
physical  barriers  of  the  hypersonic  flight  envelope  is  imposed  by  the  tremendous  heat 
loads  experienced  in  flight.  Every  undergraduate  engineering  student  is  taught  that  heat 
transfer  occurs  due  to  three  basic  mechanisms:  convection,  conduction  and  radiation.  In 
most  flight  regimes,  the  aerodynamicist  can  simply  concentrate  on  convection  and 
conduction,  completely  ignoring  the  contributions  from  radiation.  For  flows  around 
reentry  vehicles,  this  approach  is  no  longer  valid.  In  order  to  perform  analyses  of  the 
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radiating  shock  layers,  it  is  absolutely  essential  to  account  for  the  contributions  of  air 
chemistry  and  radiation  effects  to  the  overall  heat  transfer  problem.  Additionally,  severe 
thermal  and  chemical  nonequilibrium  are  known  to  occur  at  the  high  Mach  numbers 
characteristic  of  reentry  conditions,  therefore  it  is  necessary  to  account  for  these 
relaxation  processes  as  well. 

Finally,  it  is  noted  that  radiation  is  transmitted  both  towards  and  away  from  a 
body  reentering  the  atmosphere.  The  radiation  transmitted  toward  the  body  is  of  interest 
to  the  design  engineer  who  is  concerned  with  the  heat  transfer  problem  described  above. 
However,  the  radiation  transmitted  away  from  the  body  is  of  interest  not  only  to 
engineers  and  researchers  but  to  MASINT  personnel  also.  It  is  conceivable  that  a 
competitor  state  with  an  adequate  technological  and  industrial  base  could  in  the  near 
future  pursue  stealthy  hypersonic  weapon  systems  as  a  deterrent  to  future  US  systems  of 
like  construction.  Given  the  speeds  involved  in  a  hypersonic  strike,  the  time  to  detect  an 
enemy ‘s  hypersonic  weapon  system  will  be  perhaps  the  critical  link  in  the  kill  chain  for 
US  countermeasures.  Therefore,  it  is  relevant  to  the  national  security  of  the  US  to  be  able 
to  detect  these  stealthy  vehicles.  Fortunately,  while  it  may  be  possible  to  reduce  the  radar 
cross-section  of  a  vehicle  or  even  to  mask  propulsion  signatures,  it  is  impossible  to 
conceal  the  radiation  which  is  emitted  by  the  highly  energetic  gas  in  the  shock  layers 
surrounding  vehicles  moving  at  high  Mach  numbers.  The  degree  to  which  one  can 
correctly  model  these  phenomena  will  have  a  direct  bearing  upon  the  ability  of  future 
MASINT  personnel  to  correctly  identify  hostile  vehicles  en  route  to  attack  the  US,  its 
interests  and  its  deployed  forces.  It  is  in  light  of  these  considerations  especially  that  the 
method  which  is  developed  in  the  following  chapters  is  proposed. 
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II.  Background 


Computational  models  have  been  used  extensively  in  the  field  of  hypersonics 
since  the  1950s  and  1960s.  Many  of  these  early  calculations  based  on  curve-fitted 
experimental  data,  while  others  were  more  closely  related  to  first  principles.  In  either 
case,  it‘s  rather  impressive  that  these  models  provided  reasonable  enough  estimates  to 
design  the  reentry  vehicles  of  the  early  manned  space  programs  (Berman,  1983).  From 
this  starting  point,  the  evolution  of  these  numerical  models  has  naturally  followed  the 
evolution  of  the  digital  computer.  The  1980s  saw  many  impressive  calculations  from 
first  principles  for  geometrically  simplified  flowfields.  The  developments  of  the  last  two 
decades  have  brought  about  epochal  improvement  in  computational  capabilities,  and  so 
the  application  of  theory  to  computational  models  has  continued  to  advance.  These 
advancements  have  basically  followed  some  combination  of  these  trends:  1)  higher 
dimensional  flow  fields,  2)  more  accurate  physical  models  and  3)  coupling  of  physical 
phenomena  (i.e.,  flowfield,  ablation,  radiation,  material  response,  etc.). 

Today,  it  is  possible,  although  rare,  to  see  research  codes  capable  of  calculating 
nonequilibrium  flowfields  coupled  with  radiation  and  ablation  (Johnston,  2006,  2008; 
Feldick,  et  al.,  2008).  However,  even  today  the  computational  cost  of  implementing  the 
most  general  theories  in  an  aerothermodynamic  code  is  prohibitive  and  various  trade-offs 
are  made.  The  most  expensive  aspect  of  a  fully-coupled  radiation  flowfield  methodology 
is  the  calculation  of  the  spectral  data  and  radiative  transport.  This  expense  is  due  to  the 
influence  of  a  single  additional  independent  variable,  wavelength.  Radiative  transport  in 
a  participating  or  grey  medium,  such  as  a  high  temperature  gas,  depends  on  the  spectral 
properties  of  the  transport  medium  such  as  emission,  absorption  and  scattering.  (Modest, 
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2003;  Zel‘dovich  and  Raizer,  2003)  — Bud  models”  are  the  simplest  approach  to 

modeling  this  part  of  this  dependency.  These  models  are  based  upon  the  assumption  of 
equilibrium  state  populations  according  to  a  single  equilibrium  heavy-particle 
temperature.  Thus,  radiative  phenomena  for  a  given  species  are  lumped  into  a  band, 
wherein  the  magnitude  of  the  radiative  flux  is  a  direct  empirical  result  of  this  temperature 
and  the  species  number  density.  These  methods  have  the  advantage,  in  coupled  flowfield 
solutions  of  being  computationally  inexpensive.  However,  because  of  the  chemical  and 
thermal  equilibrium  assumptions  inherent  in  these  methods,  they  are  not  always 
applicable  in  hypersonic  flows  where  these  assumptions  are  violated.  (Olstad,  1971; 
Zoby,  1993)  There  have  been  some  limited  attempts  to  adjust  banded  models  for 
nonequilibrium  effects  (Greendyke  and  Hartung,  1991;  1994),  but  in  many  cases, 
radiative  phenomena  are  too  complex  to  accommodate  the  banded  models  to 
nonequilibrium  environments. 

The  radiation  observed  within  a  high-temperature  gas  exhibits  a  complex  structure 
in  terms  of  its  spectral  characteristics.  This  structure  ranges  from  the  very  coarse,  which 
spans  hundreds  or  thousands  of  nanometers,  such  as  molecular  bands,  to  the  very  fine, 
those  spanning  a  few  nanometers  or  less,  such  as  line  emission  and  absorption.  Under 
equilibrium  conditions,  the  band  models  above  account  fairly  well  for  many  of  these 
coarse  structures.  However,  the  fine  structures  are  entirely  missed.  This  lack  of  detail  is 
actually  a  rather  significant  shortcoming  since  the  finer  structures  account  for  the  bulk  of 
the  radiation  emission  and  absorption  (Herzberg,  1950).  The  hybrid  model  of  Nicolet 
(1969,  1970)  called  RAD/EQUIL  represented  an  improvement  to  the  typical  band  model. 
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This  hybrid  model  consists  of  a  mixture  of  banded  models  for  continuum  equilibrium 
radiative  processes,  and  detailed  spectral  analysis  for  certain  selected  line  emissions. 

The  most  advanced  code  for  the  calculation  of  radiative  transport  is  the  NEQAIR 
code  originally  developed  by  Whiting,  et  al.,  (1969)  as  a  spectrographic  code.  The  code 
was  later  adapted  for  thermochemical  nonequilibrium  effects  by  Park  (1985)  for  use  in 
hypersonic  flowfields.  With  this  upgraded  capability  NEQAIR,  is  now  capable  of 
calculating  the  population  of  upper  molecular  and  atomic  states  based  upon  the  heavy 
particle,  rotational,  vibrational,  and  electron  temperatures.  Following  from  the 
determination  of  the  state  populations,  the  code  then  performs  a  line-by-line  integration 
through  the  user-defined  spectral  region  under  consideration  for  the  determination  of 
local  radiative  emission  and  absorption.  NEQAIR,  while  highly  accurate,  has  a  higher 
computational  cost  associated,  relative  to  the  band  and  hybrid  models.  Johnston‘s  HARA 
code  is  another  notable  nonequilibrium  radiation  code  gaining  popularity  in  the  literature 
for  its  use  of  some  of  the  most  recent  rate  data  available  (Johnston  2006,  2008;  Feldick,  et 
al.,  2008).  The  code  to  be  used  in  this  research  effort,  SPRADIAN,  is  a  variant  of 
NEQAIR  and  was  developed  by  researchers  Fujita  and  Abe  at  JAXA  (1997).  It  has  been 
chosen  because  of  the  highly  detailed  line-by-line  method  it  uses  to  calculate  emission 
and  absorption  coefficients  and  the  ease  with  which  it  may  be  modified  to  accommodate 
the  multitemperature  thermal  model  herein  discussed. 

Despite  the  existence  of  such  a  detailed  method  as  SPRADIAN,  the  accuracy  of 
the  radiative  solution  will  only  be  as  high  as  is  allowed  by  the  accuracy  of  the  solution  for 
the  thermodynamic  state  of  the  flowfield.  Furthermore,  for  any  investigation  concerned 
with  the  radiation  resulting  from  reentry  conditions,  it  is  necessary  to  account  for  the 
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nonequilibrium  of  the  thermodynamic  state.  This  characterization  of  the  nonequilibrium 
may  be  accomplished  at  varying  levels  of  approximation.  Those  models  which  involve 
the  least  approximation  directly  simulate  the  state-to-state  transitions  among  the  internal 
energy  levels  of  the  flow  species  (Magin,  et  al.,  2008;  Josyula  2000;  Park,  1992). 

While  potentially  very  accurate,  this  class  of  models  presents  a  set  of  calculations 
that  are  far  more  expensive  to  perform  than  some  other  useful  approximations.  The 
distribution  of  energy  within  the  internal  energy  manifolds  may  also  be  approximated  by 
partitioning  the  internal  energy  modes  and  assuming  equilibrium  within  or  among  them, 
according  to  the  nature  of  the  relaxation  processes  involved.  When  the  internal  energy 
modes  are  thus  partitioned,  the  nonequilibrium  state  is  adequately  specified  with 
knowledge  of  the  species  number  densities,  Nt,  and  the  temperatures  which  are 
characteristic  of  the  thermal  nonequilibrium,  T,  Trotyi,  Tvibj  and  Te.  The  methods  of  this 
class  which  are  proposed  in  the  literature  basically  differ  according  to  the  assumptions 
made  regarding  the  nature  of  the  relaxation  processes.  For  instance,  the  popular  two- 
temperature  model  of  Park  (1985,  1992)  posits  that  the  thermodynamic  nonequilibrium 
may  be  adequately  characterized  by  utilizing  a  common  temperature  heavy  particle 
temperature  for  T  and  Tm  ,  based  on  a  heavy  particle  energy  equation  and  a  common 
electronic-vibrational  temperature  Tev  based  on  an  electron-electronic  energy  equation.  It 
has  been  proposed  that  the  nonequilibrium  may  be  modeled  with  improved  accuracy  by 
relaxing  the  second  assumption  made  under  the  two -temperature  model  (Josyula  and 
Bailey,  2003).  By  relaxing  this  assumption,  the  vibrational  temperature  is  allowed  to 
vary  by  species  and  is  no  longer  artificially  constrained  by  the  electron  temperature. 
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Thus,  this  so-called  multispecies  multitemperature  model  allows  the  internal  energy  to  be 
redistributed  in  a  more  realistic  fashion. 

In  spite  of  the  potential  benefit  of  using  the  multitemperature  model  in  a  coupled 
flow  field-radiation  solution  method,  the  two-temperature  model  pervades  the  literature 
as  the  de  facto  standard  method.  Prior  to  this  research  activity,  virtually  no  work  has 
been  performed  investigating  the  effect  of  exchanging  the  two-temperature  model  of 
thermal  nonequilibrium  for  the  multitemperature  model  in  a  coupled  flow  field-radiation 
computer  code  with  a  line-by-line  specification  of  the  radiation  transport  solution.  As 
stated  previously,  this  highly-accurate  implementation  of  the  coupled  flowfield  and 
radiation  solutions  is  important  to  the  current  methodology.  With  the  high  degree  of 
accuracy  comes  severe  computational  cost.  This  trade-off  between  accuracy  and  the 
efficiency  computation  of  solutions  is  accepted  up  front. 

The  primary  objective  of  this  research  effort  is  to  accomplish  a  loosely-coupled 
implementation  of  a  detailed  radiation  solver,  such  as  SPRADIAN,  within  a  suitable 
nonequilibrium  flowfield  solver.  Prior  to  the  early  1990s,  many  efforts  to  accomplish 
this  type  of  coupling  have  been  attempted  on  a  simplified  level  and  have  maintained  the 
equilibrium  assumptions.  A  very  notable  exception  to  this  history  has  been  the  LORAN 
code  of  Hartung  (Hartung  1991,  Chambers  1994),  which  implemented  Nicolet‘s 
RAD/EQUIL  code — utilizing  Park‘s  nonequilibrium  state  population  calculation — within 
the  LAURA  code  (Gnoffo  1990,  Cheatwood  1996).  LAURA  is  a  finite  volume  based 
method  for  nonequilibrium  hypersonic  flows  in  chemical  and  thermal  nonequilibrium 
using  finite  rate  chemical  reactions  and  Park‘s  two-temperature  model  for  thermal 
nonequilibrium  (Park,  1987).  Another  notable  exception  would  be  the  development  of 
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HARA  and  its  coupled  implementation  in  a  viscous  shock  layer  code,  again  utilizing 
Park‘s  nonequilibrium  model  (Johnston,  2006).  The  goal  of  both  the  LORAN  code  and 
the  Johnston  viscous  shock  layer  code  was  to  provide  relatively  fast  coupled  solution 
methods  for  engineering  design  and  analysis  of  hypersonic  vehicles.  To  date,  these  codes 
represent  the  ^state-of-the-art”  for  coupled  radiation  flowfield  solution  methods.  A  more 
detailed  method  would  involve  the  use  of  a  thermal  model  which  accounts  for  the 
conservation  of  vibrational  energy  on  a  species-by-species  basis,  such  as  the  method 
implemented  by  Josyula  &  Bailey  (2006).  This  computer  code,  hereafter  referred  to  as 
NH7AIR,  has  already  been  utilized  to  accomplish  simplified  uncoupled  radiative 
flowfield  calculations — both  along  the  stagnation  line  and  for  the  whole  flowfield 
(Komives,  2009  a.  and  b.;  Martin  2010). 

There  exists  only  one  detailed  radiation  analysis  code  in  practical  use  today  -  the 
NEQAIR  (SPRADIAN)  code — and  it  has  never  been  used  in  a  coupled  fashion  with  any 
flowfield  solution  method.  There  have  been  several  attempts  at  coupling  other  radiative 
solution  methods  to  flowfield  codes,  and  even  some  attempts  to  look  at  the  effects  of  both 
radiation  and  ablation  in  flowfields  (Nicolet  1970,  Sutton  1973),  but  never  with  such  a 
detailed  line-by-line  method  as  SPRADIAN.  Also,  these  research  efforts  utilized 
relatively  simple  approximations  of  the  nonequilibrium  flow  field  conditions.  The 
method  detailed  in  this  dissertation  advances  the  state-of-the-art  with  both  the  detailed 
radiation  solution  enabled  by  SPRADIAN  and  the  enhanced  nonequilibrium  flow  field 
solution  provided  by  NH7AIR. 
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This  work  reported  in  this  dissertation  was  conducted  according  to  the  following 
research  objectives.  These  research  objectives  are  intended  to  be  overarching,  with 
specific  supporting  details  deferred  to  Chapters  III,  IV  and  V. 


Research  Objective  1 :  Develop  a  computer  code  suitable  for  the  loosely- 
coupled  calculation  of  nonequilibrium  radiative  heat  transfer  within 
reentry  body  shock  layers.  The  calculation  of  emission  and  absorption 
coefficients  shall  be  performed  by  a  highly  accurate,  line-by-line  method. 
The  nonequilibrium  state  of  the  flowfield  shall  be  solved  by  via  a 
multitemperature  nonequilibrium  flow  solver  capable  of  separately 
tracking  the  vibrational  energy  of  each  diatomic  species  and  the  energy  of 
the  free  electrons. 

Research  Objective  2:  Implement  two  radiative  transport  schemes.  The 
first  shall  utilize  the  tangent  slab  assumption.  The  second  shall  be  based 
on  the  conservation  of  radiative  energy;  namely,  it  shall  be  a  finite  volume 
method  scheme. 

Research  Objective  3:  Validate  the  developed  computer  code  against  the 
benchmark  FIRE  II  flight  experiment  (Lewis  and  Scallion,  1966;  Comette, 
1966;  Cauchon,  1972). 


Chapters  III  and  IV  provide  a  summary  of  the  theory  and  methodology,  respectively, 
which  attend  the  present  research  activity.  Given  the  complexity  of  the  problem 
investigated,  a  generous  amount  of  space  has  been  devoted  to  a  discussion  of  the  theory 
and  computational  methodology.  Specific  details  regarding  the  computer  implementation 
of  the  above  theory  and  methodology  follow  in  Chapter  V.  Results  are  presented  in 
Chapter  VI,  and  the  conclusions  drawn  from  the  performed  research  are  reported  in 
Chapter  VII. 
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III.  Theory 


The  focus  of  this  dissertation  is  the  development  and  validation  of  a 
computational  method  suitable  for  calculating,  in  a  detailed  manner,  the  flow  field- 
radiation  solution  of  typical  shock  layers  in  thermochemical  nonequilibrium  about  reentry 
vehicles  of  interest.  In  this  chapter  the  basic  theory  associated  with  the  development  of 
such  a  method  is  outlined.  The  theory  is  presented  under  two  broad  headings  concerning 
those  aspects  of  the  research  activity  which  pertain  to  the  flow  field  and  radiative 
solutions,  respectively. 

Characterizing  the  Hypersonic  Environment 

In  the  strictest  sense,  the  demarcation  between  subsonic  and  supersonic  refers  to 
that  condition  wherein  the  local  free  stream  velocity  exceeds  the  local  speed  of  sound — 
that  is  to  say,  M  >  1 .  Unlike  the  easy  distinction  made  between  subsonic  and  supersonic, 
distinguishing  between  supersonic  and  hypersonic  flows  in  terms  of  a  Mach  number  is 
somewhat  arbitrary.  This  observation  should  not  be  terribly  surprising  since — ignoring 
other  effects  for  a  moment — there  are  no  sudden  qualitative  changes  in  the  behavior  of 
the  flow  relative  to  the  propagation  of  acoustical  information  within  the  domain,  as  there 
are  when  a  flow  reaches  the  speed  of  sound.  Despite  the  inherent  limitation  of  such  a 
description,  as  a  general  rule  of  thumb,  this  change  is  said  to  occur  in  air  somewhere 
around  Mach  5  (Anderson,  2006).  Yet,  it  is  perhaps  more  instructive  to  say  that  the 
hypersonic  regime  is  one  characterized  by  certain  flow  features  and  physical  phenomena 
which  become  increasingly  influential  upon  flow  behavior,  with  increasing  Mach 
number,  and  that  this  influence  is  first  appreciable  around  Mach  5. 
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Anderson  (2006)  and  Bertin  (1994)  proceed  from  this  observation  to  delineate 
some  of  these  characteristic  features  and  phenomena.  The  first  and  most  basic  such 
feature  is  that  of  the  shock  layer.  As  Mach  number  increases,  the  shockwave  moves  in 
closer  to  the  body,  and  the  air  between  the  shockwave  and  the  body  is  thus  confined  to  an 
increasingly  smaller  region.  This  thin  region  near  the  body  is  called  the  shock  layer.  It  is 
convenient  to  discuss  the  other  characteristics  as  they  become  of  significance  in  terms  of 
increasing  velocity  and  altitude.  With  increasing  velocity,  the  bow  shock  becomes 
incredibly  strong  and  the  kinetic  energy  of  the  free  steam  is  increasingly  transferred  to  the 
internal  energy  modes  of  the  gas  particles.  This  energy  transfer  leads  to  vibrational 
excitation  and  ultimately  chemical  reactions — the  dissociation  of  molecular  oxygen  and 
nitrogen  and  eventually  the  ionization  of  the  constituent  flow  species.  Also, 
nonequilibrium  chemical  and  thermal  conditions  become  significant  due  to  the  slow 
characteristic  time  scales  of  these  relaxation  processes  relative  to  the  time  scales  of  the 
flow.  The  character  of  the  flow  also  changes  with  increasing  altitude.  At  altitudes  which 
are  sufficiently  low,  the  mean  free  path  of  the  flow  is  small  enough  relative  to  the 
characteristic  length  scales  of  the  flow  that  the  continuum  approximation  may  be 
assumed.  As  altitude  is  increased,  interesting  features  begin  to  arise.  First  the  entropy 
layer  increases  in  height  and  begins  to  engulf  the  boundary  layer,  thus  introducing  a 
troublesome  vorticity  interaction  via  Crocco‘s  theorem.  Also,  the  boundary  layer  and 
shock  layer  begin  to  merge,  and  the  shock  layer  thickens.  As  altitude  increases  still 
further,  so  does  mean  free  path,  and  transition  begins  to  the  free  molecular  regime,  where 
thermal  and  velocity  slip  become  important  effects  near  the  wall.  Furthermore,  the 
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continuum  assumption  begins  to  deteriorate,  and  it  becomes  necessary  to  transition  to  an 
appropriate  particle  description  of  the  system. 

The  Knudsen  number,  Kn  =  XI L ,  is  the  ratio  of  the  mean  free  path  to  the 
characteristic  flow  length  and  is  a  parameter  used  primarily  to  distinguish  continuum 
flow  conditions  from  non-continuum  flow  conditions.  As  such,  it  provides  a  convenient 
rule  of  thumb  regarding  applicability  of  continuum  formulations  to  a  set  of  flow 
conditions.  For  Kn  <  0. 1  ^  the  flow  is  said  to  be  in  a  continuum  regime,  while  for 

Kn»  1.0  the  flow  is  said  to  be  in  the  free  molecular  regime.  For  continuum 

calculations,  the  flow  situation  may  be  calculated  from  the  Navier-Stokes  equations.  The 
non-continuum  conditions  require  using  an  appropriate  kinetic  or  particle-based, 
description  (Evans  and  Harlow,  1957;  Bird,  1994).  While  it  has  been  a  common  practice 
to  utilize  one  method  or  another  in  the  course  of  a  particular  investigation,  it  has  been 
observed  that  a  certain  hybrid  flows  may  exhibit  regions  of  transition  between  continuum 
and  non-continuum  conditions.  So-called  hybrid  solvers  have  been  proposed  and 
extensively  developed  to  provide  a  method  of  treating  these  flow  situations  (Kolobov,  et 
al.,  2006).  For  the  investigated  trajectory  points  of  the  FIRE  II  experiment,  the  Knudsen 
number  is  around  0.001  and  thus  the  continuum,  or  fluid,  description  is  applicable. 

Governing  Fluid  Equations 

Given  the  applicability  of  the  continuum  assumptions  to  the  present  investigation, 
it  is  possible  to  utilize  a  suitable  Navier-Stokes  solver.  The  Navier-Stokes  equations  in 
their  canonical  form  do  not  address  chemical  or  thermal  nonequilibrium  or  any  of  the 
relaxation  mechanisms  associated  with  these  conditions.  Therefore,  in  order  to 
accommodate  this  formalism  to  the  study  of  high-temperature  gas  flows,  in  which  these 
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features  are  important,  several  further  improvements  ‘  must  be  made  by  incorporating 
additional  source  terms  to  model  these  sources  of  nonequilibrium  .  These  improvements 
upon  the  basic  governing  equations  have  been  well  developed  by  others  and  are  presented 
below  in  the  subsections  which  follow  (Park,  1992;  Lee,  1986;  Appleton  and  Bray,  1964; 
Holt,  1965). 

Conservation  of  Species  Mass. 

In  a  reacting  multispecies  flow  it  is  necessary  to  account  for  the  effects  of  both 
species  diffusion  and  species  production  and  destruction  (Josyula  and  Bailey,  2003). 


(1) 


The  first  term  on  the  right  hand  side  represents  the  divergence  of  the  species  mass  flux 
vector,  as  illustrated  by  the  definition  of  the  diffusion  velocity. 


vj  =ut 


(2) 


And  the  second  term  on  the  right  hand  side  is  the  species  production  term.  For 
continuity,  it  is  required  that  each  of  the  two  new  terms  equal  zero  when  summed  over  all 
the  constituent  species. 

Z^=0  (3) 

!>//  =0  (4) 


It  is  also  noted  that  the  mixture  density  is  obtained  from 

p=Tjps=TjNsms-  (5) 

s  s 


The  production  term  d>s  accounts  for  the  contributions  to  species  s  by  chemical  sources 
and  sinks  in  the  flow  such  as  dissociation,  ionization,  recombination  and  attachment. 
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These  processes  are  modeled  via  the  use  of  the  rate  equation  (Vincenti  and  Kruger, 
1967).  Consider  a  gas  mixture  of  species  5  undergoing  r  elementary  chemical  reactions 


,  kku  ^Vv.x 

5=1  b'r 


(6) 


5=1 


where  v  \  r  and  v  "s  r  are  the  stoichiometric  coefficients  of  the  reaction  and  kfr  and  kb,r  are 

the  forward  and  backward  rate  constants.  Whereby,  the  contribution  of  reaction  r  to  the 
rate  of  change  of  the  concentration  of  species  5  is  given  by 


4X5] 

dt 


!]  ={v\f-y\k)kuY[\_xs-] "  +(v’v-vV)^n[Xs]  "  (10) 

Jr  s= 1  s=\ 


/ 

1 

s= 1 


In  this  way,  the  total  rate  of  change  of  the  molar  concentration  of  species  s  [Xs  ]  is  given 
by  summing  over  all  contributing  reactions  r. 


dt  s'r  s’r > 


kf,ri\[xs]  -4,nw 


s= 1 


(11) 


For  equilibrium  calculations,  the  implementation  of  the  foregoing  equations  can  be 
simplified  by  use  of  the  so-called  equilibrium  constant  from  the  law  of  mass  action, 
which  relates  the  forward  and  backward  reaction  rates  according  to 


K 


(12) 


Conservation  of  Momentum. 

The  species  conservation  of  momentum  equation  is  given  by 


|(a«;  )  ■ V  )+-£j{p,«‘KJ  +  P,»JK )  +  fj  - 


dxJ 


dxJ 


dxJ 


F'  +Fl  +Fl 

ele,s  ela,s  inela,s 


(13) 
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The  various  forces  Flm  s  arise  from  different  collision  types — namely,  elastic  or  inelastic, 

neutral  or  charged  species.  Because  of  its  importance  in  plasma,  the  electron  momentum 
equation  is  stated  explicitly  here 


r 


where  E‘  is  the  electric  field  associated  with  either  an  external  or  induced  electric  field 
and  if  is  the  effective  collision  frequency  as  given  by  Lee  (1983,  p.38).  The  electric 

field  solution  in  the  absence  of  an  external  or  induced  magnetic  field  reduces  to  a  solution 
of  the  Poisson  equation  by  which  the  space  charge  distribution  within  plasma  is  related  to 
the  electrostatic  potential.  The  numerical  solution  of  the  Poisson  equation  is  relatively 
expensive.  Various  approximate  models  may  be  used  in  order  to  evaluate  the  electric 
field  in  a  computationally  efficient  manner,  such  as  the  electron  gas  pressure  or 
ambipolar  diffusion  approximations.  In  the  current  implementation,  it  is  assumed  that 
electrons  and  ions  diffuse  in  an  ambipolar  fashion,  wherein  the  fluxes  of  electrons  and 
ions  are  assumed  to  be  equal  and  related  to  one  another  through  the  ambipolar  diffusion 
coefficient  Da 


T  =T,  =T  =  DS7N 

el  a 

The  electric  field  in  such  situations  may  be  given  by 


fVA0 

- 1 

+ 

N  J 

(15) 


(16) 
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Having  approximated  the  electric  field  in  this  way,  the  volumetric  electric  field  force  F'el 
is  given  by 


EL  =Y.N.eZ,E‘ 

s 


(17) 


Finally,  summing  up  the  contribution  from  each  species  the  total  momentum  equation  is 


given  as 


d_ 

dt 


+  —  (puluj 

dxJ  y  s 


dp 

+  -£1t  - 

dxJ 


dxJ 


T7‘  +  F'  +  F‘ 

‘  ele  T  '  ela  T  1  inela 


(18) 


Thermal  Nonequilibrium . 


The  popular  two-temperature  model  approximates  the  nonequilibrium  situation  by 
postulating  the  existence  of  a  common  heavy-particle  temperature  7  characteristic  of  the 
translational  and  rotational  energy  modes  and  a  second  temperature  Tv  characteristic  of 

the  energy  contained  within  the  combined  vibrational  energies  of  all  diatoms,  the 
electronic  energy  and  the  energy  in  the  free  electrons.  This  nonequilibrium  model  utilizes 
a  single  vibrational-electronic  energy  equation  from  which  Tv  is  calculated  (Gnoffo,  et 


al.,  1990). 


d  /  n.  d  (  :\  du‘  d 


dt 


dxJ  dxJ 


(  \  ^ v 


l  X  (Pe)s,r  +  Qele,e  +  Qinel,e  +  Qt-V  ~  Qrad  +  X  ®sD 
V  ™  J  s,r  s=mol.  s 


(19) 


The  justification  for  a  two-temperature  model  is  based  on  two  considerations:  1) 


the  rapid  energy  transfer  between  the  translation  of  free  electrons  and  the  vibrational 


motion  of  the  diatoms  and  2)  that  the  distribution  of  internal  energy  among  the  low-lying 
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electronic  states  of  the  heavy  particles  may  be  characterized  as  being  in  equilibrium  with 
the  ground  electronic  state  at  the  electronic  temperature. 

The  multispecies  multitemperature  thermal  model  referred  to  throughout  this 
dissertation  consists  in  the  following  details  regarding  the  conservation  of  energy  among 
the  various  internal  energy  modes  within  the  flow,  together  with  the  various  terms  which 
model  the  energy  exchanges  which  take  place  between  them.  The  key  feature  of  this 
multispecies  multitemperature  model  is  the  separate  tracking  of  energy  in  the  electron- 
electronic  state,  Equation  (23),  and  in  the  vibrational  energy  manifolds  of  each  diatomic 
species,  Equation  (27). 

Conservation  of  Energy. 

The  expression  for  the  conservation  of  internal  energy  for  atomic  and  diatomic 
species  is  given  by 


8  f 

j,p ■ 


-ir+u'r+e 

v2 


d  j(  l  2  ^  8qJs 

dxJ  s  l  2  s  dxJ 


d  f1  +  - —(u‘tv)  = 

dx‘ y  s  s  ' 


8xJ 


J 


dxJ 


-Ps11  Vs  +-PsUuVs 

\Z  z 

"  '  X  (Pe)s,r  +  Pele,s  +  Qele,s  +  Qinel,s  ~  Q, 


(20) 


V  N 


rad 


J  s,r 


where  the  total  species  energy  is  given  by  the  following  equations  according  to  the 
species  kind. 


Atoms: 


3  k  Tx 

e,  = - T  + 

2  m„ 


r  k' 


VmsJ 


0 


el,s 


{giJgoZ^zl 

l  +  {gxJ  g*,s)e~9eU 


(21) 


Diatoms: 


5  k 

e  = - T  + 

2  m. 


f  k^ 


\msJ 


0 


-0vJTv 


-  + 


\msJ 


0 


el,s 


(gM  /_goZK^l 

1+(gi,s/go,,)e"^s/7; 


(22) 
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The  electron-electronic  energy  conservation  equation  is  very  similar  to  the  species 


conservation  equation. 

d  ( 1 


i,p- 


—  U  +  U  Vn  + 


-I - - p  uJ  I  —  u2  +e  -I - —  + 

1  r  e  I  r\  e  I 


J 


dxJ 


2 


dxJ 


d  r 
dxj 


\p,u2v;  +  \PyuT; 

z 

M 


dxl 


(23) 


ZNesr  \Pe)s,r  +  Pele,e  +  Qele,e  +  Qinel,e  ~  Qrad,e 

JX  /  s.r 


where  the  internal  energy  of  electrons  is  given  by 

[3/2  kT,+(e^,) 

p  —  — - 

Me 


(24) 


The  electron  energy  conservation  equation  is  stated  separately  from  the  general  species 
energy  conservation  equation  here  to  illustrate  the  additional  source  terms  which  play  a 
relatively  important  role  in  the  overall  energy.  The  first  such  term  represents  the  energy 
gained  by  the  production  of  electrons. 

2X,  (25) 


The  next  term  Pele  s  models  the  work  done  by  the  electric  field  upon  the  electrons,  which 


is  also  known  as  Joule  heating. 


Pele,s=NsZsEX 


(26) 


In  lieu  of  a  more  detailed  approach,  such  as  solving  the  vibrational  master 
equations  via  a  detailed  balancing  procedure  (Park,  1992),  the  vibrational  energy 
conservation  equation  is  solved  for  each  of  the  diatomic  species  using  the  macroscopic 
nonequilibrium  vibrational  temperatures  and  the  Landau-Teller  formalism  (Landau  and 
Teller,  1936) 
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(27) 


Qt-V  ~l~  Qv-V  Qe-v  Qrad-v  +  d)sDs 


The  vibrational  energy  ev  s  represents  the  total  energy  in  the  vibrational  manifold  of  the 

diatomic  species,  which  are  modeled  as  harmonic  oscillators.  The  energy  levels  are 
assumed  to  be  populated  by  a  Boltzmann  distribution  at  the  species  vibrational 
temperatures.  This  assumption  holds  well  for  low  vibrational  states.  High  vibrational 
states  deviate  from  this  assumption,  but  the  total  energy  contained  in  these  higher  levels 
is  negligible  (Lee,  1985).  The  source  terms  Qm_v  model  the  exchange  of  energy  between 
the  various  energy  modes.  The  effect  of  vibrational  population  depletion  arising  from 
dissociation  is  accounted  for  in  the  vibration-dissociation  coupling  term  6>sDs . 

The  total  energy  conservation  equation  is  given  by 


+ 


dqj  d 
dxj  dxj 


_a_ 

dxj 


s 


Qrad 


(28) 


The  translational-rotational  temperature  is  obtained  from  equation  (29)  by  solving  for  T  . 
In  a  similar  fashion  the  species  vibrational  temperatures  are  recovered  by  solving  the 
following  expression  for  the  energy  contained  in  the  harmonic  oscillator  for  Tv  s 
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These  are  not  only  primary  thermodynamic  quantities  of  interest,  but  they  are  also  needed 
to  calculate  reaction  rates  at  the  next  time  step.  Additionally,  they  are  needed  as  input 
into  the  radiation  solver. 


Transport  Processes. 

Accounting  for  the  transport  of  mass,  momentum  and  energy  in  the  conservation 
equations  is  made  somewhat  more  of  a  challenge  for  a  high-temperature  gas  mixture.  In 
order  to  calculate  those  terms  which  account  for  the  non-convective  transport  of  these 
quantities,  it  is  first  necessary  to  compute  both  mixture  and  species  diffusivities, 
viscosities,  and  thermal  conductivities.  The  species  viscosities  are  calculated  from  the 
curve  fits  of  Blottner  (1971),  which  are  known  to  be  reasonably  accurate  up  to  10,000  K 
(Josyula  and  Bailey,  2003).  Unfortunately,  reentry  flowfields  contain  regions  that  are 
commonly  at  temperatures  well  outside  this  range.  Therefore,  these  viscosity  values  are 
used  tentatively,  and  it  is  here  noted  that  it  would  be  desirable  at  a  later  time  to 
implement  better  suited  curve  fit  data  (Gupta,  et  al.,  1987),  thus  extending  the  viscosity 
calculation  out  to  30,000  K  and  reducing  the  uncertainty  inherent  in  the  current  approach. 
The  species  thermal  conductivities  were  calculated  via  Eucken‘s  relation  (Vincenti  and 
Kruger,  1963) 


Kt=M, 


f  5 

—  C  +C 

2  v,trans,i  v,int,i 

. 


(31) 


From  the  species  transport  properties  the,  mixture  viscosity  and  thermal  conductivity 
were  computed  by  using  Wilke‘s  semi-empirical  mixing  rule  (Bird,  1960) 


XT  Mg  Pa 

[xL,«v 

P 


(32) 
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where 


-1/2 

1  + 

f  \ 

Pa. 

1/2 

1/4 

l 

l-*.J 

(33) 


Finally,  it  is  noted  that  only  ordinary  diffusion  is  accounted  for,  whereby  Ficke‘s  Law 

Piui  =  -pDyVq  (34) 


is  adequate  in  providing  an  estimate  of  the  diffusion  mass  flux  resulting  from  gradients  in 
the  species  concentrations.  This  theory  is  in  contrast  with  the  higher  order  binary 
diffusion  processes  described  by  the  Stefan-Maxwell  equations  (Cussler,  1976),  which 
account  for  the  influence  of  the  diffusion  of  other  species  upon  the  diffusion  of  the 
species  of  interest. 


Kinetic  Processes. 

The  theoretical  study  of  nonequilibrium  gasdynamics  is  in  large  part  an  effort  to 
understand  and  describe  the  nature  of  the  kinetic  processes  which  restore  a  gas  to  its 
equilibrium  condition.  Entire  volumes  could  be  (and  have  been)  filled  in  efforts  to 
catalog  the  various  models  which  have  been  advanced  to  characterize  these  processes. 
Volumes  196  and  197  of  the  AIAA  series  Progress  in  Aeronautics  (2002,  2004)  provide 
such  a  listing  of  the  most  up-to-date  information.  The  interested  reader  may  consult  these 
volumes  for  additional  details  concerning  alternate  approaches  if  desired.  Only  those 
process  models  which  have  been  adopted  in  this  solution  methodology  are  discussed 
here.  For  ease  of  discussion,  it  is  convenient  to  group  these  processes  under  one  of  the 
following  headings:  vibrational  relaxation,  chemical  reactions  and  thermal  ionization 
(Stupechenko,  1967). 
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Vibrational  Relaxation. 


The  technique  for  modeling  vibrational  relaxation  follows  from  the  theoretical 
development  presented  in  Vincenti  and  Kruger  (1967)  of  the  Landau  Teller  formalism. 
Making  use  of  this  theory,  the  rate  of  change  of  the  vibrational  energy  due  to  the 
translational-vibrational  coupling  can  be  modeled  by  a  simple  linear  ordinary  differential 
equation  of  the  form 


dE „  E*(T)  —  EV 

dt  t(T,  p ) 


(35) 


While  the  solution  to  the  above  differential  equation  is  well-known,  it  depends  on  the 
local  macroscopic  thermodynamic  state  via  the  experimentally  derived  relaxation  time 
ts(T,p).  The  species  translational -vibrational  relaxation  time  is  calculated  from  the 

Landau-Teller  interspecies  relaxation  times  zs  k  according  to 


I.*. 

_ s _ 


(36) 


The  work  of  Milikan  and  White  (1963)  furnishes  a  suitable  method  and  experimental 
data  whereby  to  approximate  the  value  of  z  k .  This  approximation  is  accomplished  by 


way  of  the  experimental  correlation 

pzsk  =  exp^(r_1/3  -0.015//1/2)-18.42j  (atm  sec)  . 


(37) 


In  this  manner,  the  translational-vibrational  energy  exchange  source  term  may  be  cast  as 

(el~ev,s) 


Qt-V  ~  Ps 


(38) 


The  electron- vibrational  energy  exchange  may  be  modeled  as  proposed  by  Lee  (1985). 
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(39) 


through  electron  impact  is  roughly  two  orders  of  magnitude  less  efficient  for  O2,  NO  and 
NO+  (Park,  1992).  However,  if  it  should  become  of  interest  to  investigate  the 
contribution  of  electron  impact  excitation  to  the  vibrational  modes  of  these  other  species 
summaries  of  cross  section  and  rate  data  may  be  obtained  from  the  work  of  Ali  (1981) 
and  Slinker  (1982). 

Chemical  Reactions. 


The  species  mass  source  terms  in  the  conservation  equations  are  intended  to 
model  the  contribution  of  the  various  chemical  reactions  involved  such  as  dissociation, 
ionization,  recombination  and  attachment.  In  order  to  calculate  these  rate  terms,  it  is  first 
necessary  to  have  a  means  whereby  to  do  so.  For  weakly  ionized  flows,  the  7-species  air 
model  is  considered  to  be  a  reasonable  approximation  of  the  significant  kinetic  processes. 
The  species  considered  in  the  seven-species  model  (which  the  flow  solver  presently  uses) 
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are  02,  0,  N2,  N,  NO,  N0+  and  e-  and  the  significant  kinetic  processes  are  expressed  in 
the  following  reaction  equations. 


N2+M  - ^  2N  +  M 

02+M  - - ^  20 +  M 

NO  +  M  - - ^  N+  O  +  M 

N2+0  - - ^  NO  +  N 

NO  +  O  - - "  02+N 

N  +  0  - - x  N0++  e 


(42) 


The  first  three  reactions  are  dissociation-recombination  and  the  fourth  and  fifth 
are  exchange  reactions.  The  reaction  rates  for  each  equation  are  calculated  according  to 
the  Arrhenius  equation  and  the  equilibrium  constant  which  have  been  extended  by  use  of 
an  effective  temperature  Ta  which  takes  the  place  of  the  usual  equilibrium  temperature 


kf(Ta)  =  Cfra^V(OdITa) 


(43) 


(44) 


The  constants  required  to  evaluate  the  forward  reaction  rates  kf  and  the  experimental 
curve-fit  for  Keq  are  taken  from  (Park,  1985),  more  recent  data  on  reaction  rates  may  also 
be  found  in  Park  (1989,  1990,  1992),  Gupta  and  Yos  (1987)  and  Bose  and  Chandler 
(1997).  The  calculation  of  the  effective  temperature  varies  according  to  the  nature  of  the 
reaction.  For  the  dissociation-recombination  reaction  Ta  =  T‘^'Tq  is  calculated  according 


to  the  empirical  relationship  proposed  by  Park  (1992)  to  model  the  vibrational- 
dissociation  coupling.  Typical  values  of  the  exponent  q  range  from  0.3  to  0.5;  0.5  is  used 
in  the  flow  solver  used  in  the  present  work.  The  reaction  rates  for  exchange  reactions 
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depend  upon  the  translational  temperatures  of  the  diatomic  molecules  only,  that  is  to  say 
Ta=T 

Thermal  Ionization. 

In  the  7-species  air  model,  there  are  two  reactions  by  which  NO  is  ionized — 
associative-ionization  and  electron-impact  ionization.  At  speeds  below  about  9  km/s  the 
associative-ionization  process  is  dominant  due  to  its  relatively  low  reaction  threshold 
(2.67  eY),  as  compared  to  that  of  electron-impact  ionization  (9.25  eY).  In  fact,  the 
associative-ionization  provides  the  seeding  electrons  which  are  responsible  for  the 
subsequent  production  of  electrons  in  the  electron  avalanche  process.  As  speed  increases, 
the  electron-impact  ionization  becomes  dominant,  and  other  species  begin  to  ionize  as 
well,  forming  N2+,  02+,  N+  and  0+.  Again,  in  the  present  investigation  only  the 
formation  of  NO+  is  considered  (Park,  1985,  1986,  1987).  The  effective  temperatures  for 

associative-ionization  and  electron- impact  ionization  are  taken  as  T  and  sTT  e , 
respectively,  although  electron-impact  ionization  is  neglected  in  the  present  study.  Only 
the  forward  reaction  rates  have  been  discussed  in  the  subsections  above.  The  reverse 
reaction  rates  may  be  calculated  by  using  the  equilibrium  constant  at  the  appropriate 
effective  temperature. 
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Radiative  Transport 

Radiation  in  gas  fields  has  long  been  studied  by  astrophysicists,  for  whom  the 
spectral  characteristics  of  the  radiation  from  stars  are  the  main  experimental  verification 
of  astrophysical  theories.  As  a  result  of  these  activities,  a  mature  body  of  theory  and  a 
wealth  of  data  exist  concerning  the  spectral  behavior  of  the  gas  species  which  must  be 
modeled  in  a  radiation-gasdynamic  solution  method  (Pai,  1963;  Hertzberg,  1950).  In 
developing  a  suitable  solution  method,  it  is  necessary  to  not  only  model  the  spectral 
behavior  of  the  participating  species  but  also  the  transport  of  radiation  within  the  solution 
domain.  Here  again,  a  mature  body  of  theory  exists  with  an  extensive  array  of  highly 
sophisticated  methods  and  techniques  available  to  model  the  transport  of  radiation  in  a 
participating  media,  such  as  the  high-temperature  air  surrounding  a  reentry  vehicle 
(Modest,  2003).  Therefore,  the  challenge  in  developing  a  suitable  solution  method  lies 
not  in  an  inadequate  theoretical  basis  for  the  proposed  models,  rather  in  the  fact  that  the 
computation  of  the  most  general  transport  models  is  prohibitively  expensive.  As  such, 
varying  degrees  of  approximation  are  accepted  in  the  solution  method — typically  with 
regard  to  either  or  both  the  determination  of  the  spectral  behavior  of  the  media  or  the 
solution  of  transport  equation  within  the  media.  For  instance,  in  years  past  step  radiation 
models  were  quite  standard  in  modeling  the  spectral  behavior  of  the  media,  while  in 
modeling  the  radiation  transport,  it  has  been  a  common  practice,  even  to  the  present,  to 
use  the  tangent-slab  approximation.  In  order  to  address  some  of  these  issues,  the 
following  subsections  will  address  in  turn  the  various  sources  of  emission  and  absorption 
within  a  high-temperature  gas,  the  transport  of  the  resulting  radiation,  and,  finally,  the 
coupling  of  the  radiative  energy  into  the  governing  gasdynamic  equation  set. 
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Emission  and  Absorption  Mechanisms. 


Calculating  the  solution  of  the  radiative  transport  equation  is  a  relatively 
straightforward  matter  in  cold  air,  vis-a-vis  high-temperature  air.  The  reason  that  the 
latter  is  so  much  more  difficult  is  that  media  such  as  high-temperature  air  are  said  to 
participate  ‘  in  the  transport  calculation  by  emitting  and  absorbing  radiation  within  the 
solution  domain;  this  transport  is  in  addition  to  the  radiative  transport  which  is  a  result  of 
the  radiative  boundary  conditions.  Therefore,  where  radiative  transport  in  high- 
temperature  gases  is  concerned,  two  major  tasks  are  presented  to  be  accomplished.  The 
first  task  concerns  the  determination  of  the  so-called  emission  and  absorption  coefficients 
within  the  gas  volume.  (The  second  task  is  the  subsequent  solving  of  the  radiative 
transport  equation.  This  discussion  is  deferred  to  the  next  section.)  These  coefficients 
are  calculated  through  a  combination  of  empirical  and  theoretical  considerations  which 
pertain  to  the  mechanisms  by  which  these  phenomena  occur — namely,  the  various 
electronic  transitions  which  are  possible  in  atomic  and  diatomic  systems.  There  are  three 
basic  classes  of  transitions  which  are  generally  discussed:  1)  free-free  (Bremsstrahlung), 

2)  bound-free/free-bound  (photoionization/-recombination),  and  3)  bound-bound  (line 
and  band  spectra).  The  methods  available  for  calculating  the  spectra  which  result  from 
the  above  transition  types  are  be  discussed  in  the  subsections  which  follow.  In  this 
presentation,  discussion  is  restricted  to  that  theory  which  is  necessary  to  establish  a 
common  understanding  of  the  basic  principles  implemented  in  SPRADIAN,  the  radiation 
solver  utilized  in  this  investigation  (Fujita  and  Abe,  1997).  Many  excellent  supporting 
texts  also  exist  and  have  been  consulted  where  they  facilitate  the  discussion  (Zeldovich 
and  Raizer,  1967;  Hertzberg,  1950;  Penner  and  Olfe,  1968;  and  Park  1992). 
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Free-Free  Transitions. 


In  high-temperature  air,  such  as  is  observed  at  reentry  speeds  (T  >  15,000if), 
atomic  species  are  much  more  prevalent  than  the  diatomic  species,  which  have  at  this 
point  mostly  dissociated.  Additionally,  the  reactions  precipitated  by  electron-impact  are 
well  underway,  whereby  electrons  collide  with  heavy  particles  such  as  ions  and  neutral 
particles.  In  these  inelastic  collisions,  the  electron  may  excite  the  heavy  particle,  possibly 
causing  the  heavy  particle  to  ionize  further,  or  alternately,  to  recombine  with  ions  to 
reduce  their  degree  of  ionization,  or  even  attach  to  neutrals  to  form  negative  ions.  These 
types  of  collision  are  not  considered  for  the  moment,  and  instead  the  reader4  s  attention  is 
directed  to  the  inelastic  collisions  between  electrons  and  heavy  particles,  particularly 
ions,  wherein  an  electron  interacts  with  an  ion  but  does  not  result  in  either  ionization  or 
recombination,  rather  in  the  emission  or  absorption  of  a  photon 

Az  +el+hv^Az  +  eu  (45) 

where  Az  is  a  heavy  particle  with  z-valence,  <?/  and  e„  are  free  electrons  in  lower  and 
upper  kinetic  energy  states  respectively,  and  hv  is  the  photon.  These  types  of  collision 
are  associated  with  the  so-call  free-free  transitions,  and  they  may  result  in  the  emission 
or  absorption  of  a  photon  with  energy  E  =  hv  due  to  the  acceleration,  W ,  experienced 
by  the  electron  in  the  field  of  the  heavy  particle  and  corresponding  to  the  resulting 
energy  exchange 

AE  =  —  ^  f  w 2vdv  wv  =  —  f  w  (t)e~i2m/‘dt  (46) 

2  c3 1  v  In  W 

In  the  literature,  the  radiation  resulting  from  this  mechanism  is  referred  to  as 
bremsstrahlung,  that  is  -braking  radiation”  (Zeldovich  and  Raizer,  1967). 
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The  spectral  coefficients  may  be  conceptualized  by  utilizing  the  idea  of  a  collision 
cross-section  cr^1  (A,Tel) .  In  the  expressions  below,  for  an  ionized  gas,  with  an 

equilibrium  electron  velocity  distribution,  the  emission  and  absorption  coefficients  can 
be  expressed  in  terms  of  electron  and  ion  number  density,  wavelength,  and  electron 
temperature 


(47) 


(48) 


However,  with  the  particular  values  of  the  cross-sections  unknown  for  the  moment,  it  is 
necessary  to  obtain  them  by  some  means. 

Consider  the  Planck  equation,  which  specifies  the  spectral  radiation  density 
distribution  for  equilibrium  radiation 


From  here,  an  expression  for  the  absorption  coefficient,  av ,  of  the  bremsstrahlung 


absorption  per  ion  per  electron  is  sought.  In  this  sort  of  collision  ions  are  considered 
stationary  and  electrons  moving  with  velocity  V  .  Then  the  absorption  of  radiation  at 
equilibrium  conditions,  in  the  frequency  interval  v  to  v+dv,  per  unit  time  per  unit 
volume  by  electrons  with  velocity  in  the  range  V  to  V+dV  is  given  by 

KN.U./V  ■  cf  (v)dv  •  a,  (l  -  ) .  (50) 
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where  the  factor  (l-e  ,1v/kTj  accounts  for  induced  emission.  The  radiation  emitted  under 
the  same  conditions  is 


N+Nef(v')dw'dqv(v').  (51) 

Noting  that  vdv  =  v'dv'  ,  the  general  relation  between  av  and  the  differential  radiation 
cross  section,  d<rv  is  given 


a 


V 


cV'2  dajy) 
8;rvV  dv 


(52) 


Finally,  noting  that  dqv  =  hvdcrv  and  by  utilizing  the  approximation  for  dqv  given  by 
Landau  and  Lifshitz  (1962) 


dqv  = 


32;r2  Z2e6 
3^3  m2c3v2  V 


(53) 


an  expression  for  the  unit  absorption  coefficient  av  in  terms  of  the  degree  of  ionization 


Z,  the  electron  velocity  v  and  the  frequency  v  is  obtained 

_  An  ZV 
3\f3  hcm2vvz 


(54) 


In  order  to  obtain  the  absorption  coefficient,  kv,  of  the  bremsstrahlung  radiation  at  the 
electron  temperature,  the  above  expression  may  be  multiplied  by  N+  and  Ne  and  averaged 
over  the  electron  velocity  by  use  of  the  Maxwell  velocity  distribution  function. 
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(55) 


Again,  the  emission  and  absorption  coefficients  may  be  related  through  Kirchhoff  s  law. 
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-hv/kT 


(56) 
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The  above  derivations  used  to  develop  expressions  from  the  emission  and 
absorption  coefficients  utilize  classical  mechanics  in  their  approach  (Griem,  1964). 
Predictions  based  on  the  classical  mechanical  approach  differ  from  the  quantum 
mechanical  calculations  by  a  factor  referred  to  as  the  Gaunt  factor  (ZeTdovich  and 
Raizer,  1967) 


(57) 


The  Gaunt  factors  used  by  SPRADIAN  are  taken  from  Peach  (1970). 
Bound-Free  Transitions. 


Now,  the  reader4  s  attention  is  turned  to  the  ionization-recombination  collisions  as 
developed  in  Zeldovich  and  Raizer  (1967).  In  these  collisions,  energy  must  be  absorbed 
or  released  by  the  electron-ion  system  as  the  electron  is  captured  or,  alternately,  freed. 
This  capturing  or  freeing  of  an  electron  may  occur  by  transferring  energy  to  or  from  the 
internal  structure  of  the  ion,  a  third  body  or  a  photon 

A'f'  +hv  —  A~u  +  e  (58) 

9 

where  Az  is  a  heavy  particle  with  z-valence,  e  is  the  free  electron  and  hv  is  energy  of  the 
photon.  Here  the  transfer  of  energy  by  a  photon  is  considered.  Let  us  begin  by 
considering  the  energy  levels  En  of  a  hydrogenic  atom  which  are  given  in  terms  of  the 
principle  quantum  number  n  which  ranges  from  1  to  co  at  the  ground  and  free  states, 
respectively, 


T  7 

T7  —  ZMZl 

^ n  2 

n 


and 
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(60) 


O  2  4 

T  _  e  _  2 k  me 
h~2^~  h2 

Here  1  denotes  the  ionization  potential,  whereby  the  binding  energy  of  an  electron  in  the 
n-th  quantum  state,  En  =  \En\  =  I  /  n2  is  obtained.  The  energy  of  the  emitted  or  absorbed 


photon  then  is  given  by 


,  „  .  „  .  mV  IrrZ 

hv  =  E+  E  = - h  -  - 


(61) 


Again,  it  is  possible  to  calculate  the  emission  and  absorption  coefficients  by 
beginning  with  the  principle  of  detailed  balancing.  The  number  of  electrons  captured  in 
the  photo-ionization  process  with  electron  speeds  in  the  range  v  to  v+dv  and  being 
captured  into  the  n-th  energy  level  of  the  ion  per  unit  volume  per  unit  time 

N+NJ(v)dv-v-t Tcn.  (62) 

where  N+  and  Ne  are  the  number  densities  of  ions  and  electrons  and  ocn  is  the  capture 
cross  section  into  the  n-th  level.  This  process  results  in  the  emission  of  photons  in  the 
frequency  range  v  to  v+dv. 

Accounting  for  induced  emission  as  before,  the  photoionization  process  from  the 
n-th  quantum  level  by  photons  of  frequency  v  in  the  range  v  to  v+dv  per  unit  volume  per 
unit  time  is  given  by 


N,  ^  f(v)dv  •  c  •  < 7„„  (l  -  ekm ) .  (63) 

where  <rvn  is  absorption  cross  section  for  a  photon  hv  into  the  n-th  state,  N„  is  the 

number  of  such  atoms  per  unit  volume.  Assuming  complete  thermodynamic  equilibrium, 
it  is  permissible  to  substitute  the  Maxwellian  velocity  distribution  for  f(\)  and  thus  arrive 
at  an  expression  for  N„ 
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(64) 


N„=N,—  e  kT  =N^ekT' 


For  a  hydrogenic  atom  the  degeneracy  of  the  n-th  energy  level  may  be  given  by 
gn  =  In2 .  In  the  expression  above  Ni  represents  the  number  of  the  atoms  in  the  ground 
state  with  degeneracy  gi  =  2.  The  excitation  energy  En  of  a  given  state  is  given  by 
En  -El  =IHZ2  (l-l/n2)  =  l(l-l/n2^J .  In  order  to  relate  the  cross  sections  for 
photoionization  and  recombination,  Saha‘s  equation  is  also  needed 


N+N  JljumkT 

=  2  - t — 

N  \  h2 


N 

where  N  -z — 1  and  where  z  and  z+  =1  are  the  partition  functions  of  the  atom  and  ion, 
Si 

respectively.  Finally,  equating  the  rates  given  by  equations  (62)  &  (63)  and  performing 
some  algebraic  substitutions,  the  following  expression  relating  the  photoionization  and 
radiative  capture  cross  sections  is  obtained  as 


z  f  mvc 

^vn  = -  —j -  acn  • 

g,\hv  ) 


The  emission  cross  section  is  given  by 


128;r4  Z4g10 

3-J3  mcih4v2v 


and,  finally,  the  absorption  cross  section  is 


64^-4  ewmZA  n  ( v 
cr  =  — — ^-r  oc  —  -2 


3^3  h6cv\5  Z2[v 


With  the  cross  section  thus  obtained,  the  absorption  coefficient  follows  according  to  the 


following  expression,  with  the  emission  coefficient  again  obtained  via  Kirchhoff  s  law 
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(69) 


oo 

v=2X°-. 

n* 

At  each  wavelength,  the  photon  energy  to  a  specific  bound-free  transition  is  calculated  by 
subtracting  the  threshold  energy  from  the  energy  at  a  particular  wavelength.  The  Gaunt 
factor  is  obtained  by  interpolation  from  the  look  up  table  as  before  (Zeldovich  and 
Raizer,  1967). 

Atomic  Bound-Bound  Transitions. 

The  line  spectra  of  atomic  systems  are  the  result  of  electronic  transitions,  whereby 
a  photon  may  be  either  emitted  or  absorbed. 

A,  +hv^Au  (70) 

Emission  of  radiation  from  an  atomic  system  occurs  in  both  a  spontaneous  and  an 
induced  manner.  Spontaneous  emission  of  a  photon,  with  wavelength  Xu! ,  occurs  due  to 

the  random  transition  of  an  electron  from  the  upper  to  lower  electronic  energy  levels  u 
and  /,  respectively.  This  probability  of  transition  is  quantified  through  the  Einstein 
coefficient  Aul  and  related  to  the  emission  coefficient  through 


The  contribution  of  the  induced  emission  is  taken  into  consideration  with  the  total 


absorption,  since  both  occur  in  proportion  to  the  incident  radiation.  The  Einstein  B 
coefficient  is  used  to  quantify  these  processes  and  is  related  to  the  net  absorption 
coefficient  through 

hi 

=— (”A-"A)®J  (72) 
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It  is  possible  to  recast  the  net  absorption  coefficient  in  terms  of  the  Einstein  A  coefficient 
through  the  relations 

(73) 

\l^ul 


Su^ul  Sl^lu  • 


The  resulting  expression  for  the  net  absorption  coefficient  is  thus 


Kx  =  ni 


AiKiSu  i,  nuSi 
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The  line  spectra  emitted  due  to  the  bound-bound  transitions  within  an  atomic 
system  are  not  discretely  distributed  over  the  wavelength  domain,  but  in  fact  have  a 
distribution  about  the  wavelength  Aul  denoted  by  above,  which  is  the  Voigt  line 

profile.  This  profile  is  the  result  of  a  phenomenon  called  line  broadening.  Three  major 
classes  of  line  broadening  are  customarily  considered:  natural  broadening,  Doppler 
broadening  and  pressure  broadening.  Natural  line  broadening  is  the  result  of  the 
uncertainty  principle  and  may  be  quantified  relatively  easily  from  an  application  of 
classical  mechanics  to  a  damped  oscillator,  whereby  it  is  found  that  for  radiation  damping 
the  half  width  of  the  line  is  independent  of  photon  wavelength  as  is  given  by 

A^=^  =  ^r  =  1.18xlO-4A.  (76) 

a>l  3mc2 

Doppler  broadening  occurs  due  to  random  thermal  velocity  of  the  radiating  atoms  and 
results  in  a  Gaussian  distribution  with  a  half  width  given  by 


A,  _  2kT\n2)  . 

Me2  V 
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Pressure  broadening  is  due  to  collisions  with  other  particles.  The  calculation  of  the 
associated  broadening  parameters  depends  on  the  type  of  collision  partner  involved. 
Stark  broadening  is  due  to  collisions  with  charged  partners.  The  resulting  collision  width 
is  calculated  according  to  the  following  empirical  relations.  The  results  of  Griem  (1964) 
are  used,  when  available,  to  quantify  the  full-half  Stark 
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(78) 


Where  data  is  not  available  from  Griem,  the  values  of  y  and  n  are  evaluated  according  to 
the  method  proposed  by  Arnold  and  colleagues  (1980) 


y  =  0.42 


r  i  y 

Broadening  also  occurs  due  to  non-resonant  collisions 
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In  the  case  of  resonant  collisions,  the  citation  within  the  source  code  comments 
matches  the  implemented  code.  From  Traving, 
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Noting  that 


co0  =2  K 


A 


(82) 
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and  taking  the  oscillator  strength  from  Ricther  (1968) 


£=1.50 


V  Si  J 
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(83) 


the  desired  line  broadening  term  is  thus  obtain 
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The  broadening  due  to  each  of  the  foregoing  mechanisms  has  been  characterized  by  a 
single  parameter,  namely,  the  full  width  at  half  maximum  (FWHM)  associated  with  each. 
However,  it  is,  in  general,  a  difficult  matter  to  describe  the  resulting  line  profde  which 
results  from  the  combination  of  these  various  widths.  Therefore,  it  is  done  in  an 
approximate  way  according  to  the  method  proposed  by  Olivero  and  Longbothom  (1977), 
whereby 
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Diatomic  Bound-Bound  Transitions . 


The  band  spectra  of  diatomic  systems  are  significantly  more  complex  to  calculate 
than  the  line  spectra  of  atomic  systems.  This  increased  level  of  complexity  arises  from 
the  fact  that  while  the  internal  energy  of  a  molecule  likewise  depends  on  its  electronic 
state,  there  is  also  a  complex  dependence  on  the  excitation  of  the  additional  vibration  and 
rotational  energy  modes  and  the  various  transitions  possible  among  them.  This  additional 
complexity  enters  through  the  calculation  of  the  Einstein  A  coefficient 


1  26;r4  1  \R,V,  SJtJ. 
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via  the  electronic  transition  strength  Rv,  J  and  the  Honl-London  factor,  Srr ,  which 


quantify  the  probability  of  vibrational  and  rotational  transitions,  respectively.  The 
electronic  transitional  dipole  moment  is  found  from  the  inner  product  of  the  vibrational 
and  electronic  wavefunctions,  and  x¥e ,  with  the  electronic  transition  dipole  moment 
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where  the  inner  product  of  the  electron  wavefunction  with  the  electron  transition  dipole 
matrix  has  been  combined  in  the  electron  dipole  moment,  Re(r ) 

R,(r)  =  (88) 


This  expression  can  be  approximated  through  the  Franck-Condon  principle,  where  the 
weak  dependence  of  the  dipole  moment  Re(r) on  the  intemuclear  displacement,  r,  is 
approximated  though  the  introduction  of  the  v‘-v“  transition  centroid  defined  as 
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(89) 


J  vJr 
J %^dr 


Since  the  vibrational  and  rotational  wavefunctions  are  independent  of  the  electronic 
coordinates  (Zeldovich  and  Raizer,  p.  317,  1967),  it  is  permissible  to  pull  /?,(/;v„)  out  of 
the  integral  such  that  the  transition  strength  is  approximated  as 


(90) 


where  the  Franck-Condon  factor  is  defined  by 


q 


v'v" 


(91) 


SPRADIAN  utilizes  a  lookup  table  to  evaluate  the  electronic  transition  strength 


the  values  of  the  table  were  computed  a  priori  according  to  the  method  outlined  by  Fujita 
and  Abe  (1997).  In  this  method,  Re(rv,v„)  is  assumed  to  be  given  experimentally  or 

through  some  suitable  quantum  mechanical  calculation,  while  the  Frank-Condon  factor  is 
calculated  as  follows. 

It  is  clear  from  equation  (91)  that  the  chief  difficulty  in  computing  qv,v„  lies  in 

obtaining  the  inner  product  of  the  vibrational  wavefunction  in  the  r  coordinate  system. 
This  wavefunction  must  be  obtained  via  a  suitable  numerical  solution  (Cooley,  1961)  to 
the  radial  Schrodinger  equation, 


h  g2vFv 
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(92) 


The  method  of  Rydberg,  Klein  and  Rees  (1931,  1932,  and  1947,  respectively)  is  utilized 
to  quantify  the  potential  energy  function  in  terms  of  the  vibrational  energy  G(v)  at  a 
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particular  vibrational  level  v  and  on  the  interval  \rinner  (v),  router  (v)] .  The  turning  points  of 
the  vibrational  motion  are  calculated  according  to  the  following  relations 


/• 

inner 


w=Jfg+/w!-/w. 


* . (v)=JtS+/(”)!  +/w. 


s(v) 

z’/  \  4h  } 

nv)  =  —  nrr  J 


dv' 


In^ljuc  Jl/2  ^G(v)-G(v')  ’ 


(93  a.) 


(93  b.) 


(93  c.) 


and 


g(v)  = 


Insjljuc  r  B„,dv' 


I 


(93  d.) 


#  -i/2  ^/G(v)-G(v') 

where  5  ,  is  the  rotational  constant  for  the  v‘  vibrational  level.  For  a  singlet  state,  the 
vibrational  energy  G(v)  may  be  given  by  the  Dunham  expansion 

z  ^  \  ( i  V  f  i  Y  (  i  V 

G(V)  =  ®,  -  +  V  -(0eXe  -  +  V  +  -  +  V  -  +  V  (94) 

J  J  J  J 

where  the  vibrational  constants,  coe,  coexe,  coeye  and  coeze,  are  taken  from  complied 

spectroscopic  data  such  as  that  found  in  (Hertzberg,  1950;  Jaffe  1987).  Interpolation  and 
extrapolation  is  based  on  a  Morse-type  function,  as  in  (Gilmore,  1992).  First,  the  inverse 
Morse  function  is  defined  by 

L{r)  =  -p{r-re),  (95) 

where  re  is  the  equilibrium  bond  distance  and  [5  is  given  as 


Yin1  c1  fj, 


(96) 
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The  inverse  Morse  function  has  the  follow  relationship  to  the  potential  energy 


U(r)  =  D, 


(97) 


where  De  is  the  quantum  well  depth  for  the  given  electronic  state.  With  the  potential 
energy  determined  as  above,  it  is  introduced  into  the  radial  Schrddinger  equation  in  order 

to  solve  for  the  vibrational  wavefunctions  y/v  and  finally  the  transition  strength  |i?vV|2 . 


The  Honl-London  factor  Srj,  quantifies  the  relative  probability  of  transition 

between  rotational  levels  J'  and  Only  singlet  bands  (i.e.,5"  =  5"'  =  0)  are  considered 
in  the  present  study.  A  compilation  of  Honl-London  factors  is  available  for  many  types 
of  transitions  (Shadee,  1964),  and  those  for  the  singlet  bands  are  listed  below  in  Table  1. 
The  current  version  of  SPRADIAN  includes  upper  to  lower  transitions  for  the  E-E,  n-n, 
n-E  and  E-n  electronic  configurations.  Line  shape  and  line  widths  are  calculated  as  in 
the  previous  section  for  the  J=0  line  of  each  band.  This  line  shape  is  stored  in  an  array 
and  used  for  the  other  lines  within  the  same  band  (i.e.,  those  transitions  which  share  the 
same  v‘  and  v“). 


Table  1.  Honl-London  Factors  for  Singlet  Band  Spectra 


A'-A" 

J'-J”=+ 1 

(P-Branch) 

J'-J”  =  0 

(Q-Branch) 

J'-J"  =  - 1 

(R-Branch) 

+  1 

(/'+1  — A')(/'+2  — A')  (  J '+  A ')  (  J '+ 1  -  A ')  (2  J '+ 1) 

(J'+A')(J'-1  +  A’) 

2(J'+1) 

2  J' 

o 

(/'+1  + A ')(/'+ 1-A') 

(2J  '+l)A'2 

(  J'+ A')(/'— A') 

(j'+i) 

J'(J'+1) 

J' 

-1 

(/'+1  + A ')(/'+ 2  + A')  ( J 

'-A')(J'+1  +  A')(2J,+  1) 

(J'-A')(J'-1-A') 

2(  J'+l) 

2J'(  J'+l) 

2  J' 
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Nonequilibrium  State  Populations  in  Atomic  Systems. 

In  the  foregoing  discussion  of  the  various  radiation  mechanisms,  it  was  assumed 
that  the  electronic  states  of  the  atomic  state  were  populated  according  to  the  equilibrium 
distribution.  This  assumption  was  only  a  convenience  of  presentation,  since  ultimately 
one  is  interested  in  the  effects  of  the  nonequilibrium  distributions  upon  the  respective 
radiation  mechanism  above  and  the  transport  of  radiation  within  the  flowfield.  In  a 
general  sense,  the  nonequilibrium  state  populations  must  be  calculated  by  time- 
integration  of  the  following  set  of  differential  equations 


r)N  171  m 

—L  =  Y,K(i,j)XlN,  +  K(c,i)NX  +  '£A(ij)NI  +A(c,i)N,N,  - 

Ot  J=i  /= i 

7  7  (98) 

m  m 

Y,K(iJ)N,N,-K(i,c)N,N,-Y,A(iJ)Nl-A(i,c)Nl. 

M  j=i 

However,  such  a  calculation  is  rather  impractical,  not  to  mention  unnecessary  at 
conditions  of  practical  interest.  Instead,  the  calculation  of  the  non-Boltzmann 
distributions  is  performed  according  to  the  Quasi-Steady-State  (QSS)  approximation,  as 
described  by  Park  (1992).  The  QSS  approximation  is  based  on  the  fundamental 
assumption  that  the  respective  sums  of  the  ingoing  and  outgoing  rates  of  transition 
between  electronic  states  are  each  much  greater  than  the  time  rate  of  change  of  the  given 
state  population 


dNt 

~dt~ 


« 


a  n  in 

'^K(i,j)NjN,+K(c,i)NX  +  Y.A(Uj)NJ+A(c,i)NtN, 


i= 1 


7=1 


(99) 


and 


dN( 

dt 


«’ZK(iJ)N,N,+K(i,c)NlN,  +  YdA(ij)N,+A(i,c)Nl. 

7=1  7=1 


(100) 
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Consequently,  one  may  set  the  LHS  of  equations  (99)  and  (100)  to  zero  and  proceed  to 
develop  a  suitable  method  for  obtaining  a  solution  to  the  set  of  m  algebraic  equations 


YJ^i,j)NjNe+K(c,i)N+N2e+Y,A(i’j)Nj+A{c’i)N+Ne 


i= i 


j= i 


(101) 


YjK(i,j)NiNe-K(i,c)NiNe-YjA(i,j)Ni-A(i,c)Ni=0. 

M  7=1 


Before  proceeding  with  the  development  of  the  master  equation,  the  author  shall 
pause  to  introduce  a  couple  of  definitions  related  to  the  ratio  of  number  densities  for  a 
particular  electronic  state  i  and  the  total  number  density  of  a  particular  species  a: 

Pt=Nl/NlE  (102a.) 


X  =  NJNaE. 


(102  b.) 


In  these  expressions,  E  denotes  the  hypothetical  equilibrium  values  of  these  number 
densities  as  given  by 


X,B  = 


g,N+Nt 


(  1.2  f  T-T  7-.  A 


2  Q+ 


h 


K2nmekTe  y 


exp 


v  *Te  / 


(103  a) 


where 


rn 

Ke=2n,e. 


i= 1 


(103  b) 


Now,  by  substitution  of  these  ratios  into  equation  (100),  the  desired  final  form  is  obtained 


I 

j= 1 


K(ij)+N’E  A(jJ) 


",e  N. 


pj+K(i,c)  +  A(c,i ) 


N 


N; 


iE 


YJK{ij)  +  K(i,c)  + 

j= 1 


N 


(104) 


Pi 
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One  must  also  take  care  that  any  solution  to  the  foregoing  equation  does  not  violate  the 
basic  conservation  relation 


NiE  NaE 

—  p'=XlT 


(105) 


Unfortunately,  this  last  expression  creates  a  system  of  m+1  equations.  With  the  help  of 
Park‘s  intuition  that  the  QSS  approximation  is  least  likely  to  be  satisfied  by  the  ground 
state,  thus  the  foregoing  derivation  results  in  a  convenient  linear  system  in  p  and/ 

Mp  =  C  +  T>x-  (106) 


where  the  matrix  M  and  the  vectors  C  and  D  are  strictly  functions  of  the  electron 
temperature  and  number  density.  For  clarity,  the  form  of  MC,  =  0 ,  C  and  D  is  illustrated 
below. 


-First  row: 


M(\,j)  =  NtEIN, 

(107  a.) 

Ct  =0 

(107  b.) 

Dt=NaC/N 

(107  c.) 

-Diagonal  elements  of  M  matrix: 


M(i,i)  =  YJK(i,j)  +  K(i,c)  + 

7=1 


A{i,j)  +  A(i,c) 


N„ 


(107  d.) 


-Off-diagonal  elements  of  M  matrix: 


^je  Mjj_  Q 

NiE  Ne 


(107  e.) 


-Vector  C,  i  ^1: 
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(107  f.) 


Cj=K(j,c)  +  A(c,j)^. 

-Vector  D,  i  ^  1: 

Dj=  0.  (107  g.) 

Nonequilibrium  State  Populations  in  Diatomic  Systems. 

The  problem  of  calculating  the  nonequilibrium  electronic  state  populations  in 
diatomic  systems  is  very  similar  to  the  problem  just  covered  for  atomic  systems.  Again 
the  most  general  solutions  would  require  the  time  integration  of  the  master  equation 

r)N  m 

-j-  =  N^K(v\v)N,.+K(c,v)NANBN,- 

v=0  .  (108) 

m 

NXK(v',v)N,-K{v,c)N,N , 

v=0 

However,  one  may  again  apply  the  QSS  assumption,  thereby  approximating  the  problem 
with  a  set  of  algebraic  equations.  To  begin,  consider  the  hypothetical  equilibrium 
number  density  of  the  electronic  state  i  of  the  diatom  produced  by  the  reaction  of  atomi 
with  atom2 


(109) 

(110) 


where  the  dissociation  limit  is  denoted  by  D  and  the  electron  temperature  by  Te  The 
QSS  master  equation  has  the  same  form  as  for  atomic  systems, 

M  =  C  +  Dj  (111) 


although  the  definitions  of  the  matrix  M  and  the  vectors  C  and  D  have  changed  as 
detailed  below 
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-First  row: 


M(i,j)  = 


2>, 


(112  a. 


C,  =0 


(112  b. 


-Diagonal  elements  of  M  matrix: 


(112  c.) 


M(i,i)  =  -K,(i,c)^-Kl(i,c)-'£  K,(i,j)^+K,(i,j) 


A(Uj) 


(112  d.) 


-Off-diagonal  elements  of  M  matrix: 


nE,eU)  nE,h{ 0  nh 


nh  nE,h  (l) 


(112  e. 


-Vector  C,  i  ^  1: 


>  =  -Ke(i,c)^Y-Kh(i,c). 


nE,h  (0 


(112  f.) 


-Vector  D,  i  ^  1: 


D(j)  =  0. 


(112  g.) 


In  the  expressions  above,  the  heavy-particle  excitation  rate  coefficients  are  calculated 
from  empirical  curve-fit  data  expressed  in  the  Arrhenius  form 


^  t  \n  It 
K,  (i,j)  -  A  — ^  exp 
kK  ’  6000 J  T 


(113) 


where  the  parameters  A,  n  and  Td  are  read  into  the  computer  program.  The  bound-bound 


electron  impact  excitation  rate  coefficients  are  calculated  numerically  as  proposed  by 
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Park.  The  average  excitation  rate  corresponding  to  the  transition  between  the  electronic 
states  e  and  e  ‘  is  expressed  as 


K{e,e')  = 


Zvexp[-^JZ^(2J+1)exp(-^ 

y  exp  -Q-  y  (2J  +  l)expf-— 
kT  )A-‘A  ’  I  kT 


£“pf-§]^(2y+1)ex[fff 


(114) 


where  C,  represents  the  electron-impact  cross  section  for  the  diatomic  species  under 
consideration.  This  expression  is  in  effect  an  average  of  the  electron-impact  transition 
coefficient  between  states  (e,v,J)  to(V,v over  the  quantum  numbers  (v,J)  as 
given  by 


K(e,v,J;e',v',J')  - 


2  nkT 


|  crexp 


(115) 


The  average  is  taken  into  account  according  to  three  considerations.  First,  the  rate 
coefficient  is  weighted  according  to  the  _mul1iplicity‘  of  each  vibrational  and  rotational 
level  within  the  initial  electronic  level.  This  weighting  is  accomplished  through  the  sums 


ZvexP 


(116) 


2^(2J  +  l)exp 


(117) 


which  are  normalize  by  the  denominator,  the  product  of  the  sums  of  these  initial  state 
multiplicities.  This  result  is  easily  confirmed  by  the  definition  of  the  partition  function 
for  the  vibrational  and  rotational  modes.  Second,  the  rate  coefficient  is  weighted 
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according  to  probability  of  a  v-to-v‘  transition,  through  the  Franck-Condon  factor, 
q(v,v') .  Finally,  the  rate  coefficient  is  weighted  according  to  the  degeneracy  of  the  final 
rotational  state,  given  by  the  quantum  number  J‘. 

Radiative  Transport  Equation 

Having  discussed  the  theory  associated  with  the  first  major  task  in  calculating  the 
radiation  field,  the  remainder  of  the  discussion  is  directed  towards  obtaining  a  solution  to 
the  radiative  transport  equation.  The  radiative  transport  equation,  in  a  non-scattering 
medium,  is  given  in  differential  form  along  a  single  ray  as 

— —  =  £  —  k'  I  (118) 

dx 

If  one  considers  a  region  where  the  emission  and  absorption  coefficients  are  uniform, 
such  as  at  equilibrium,  the  radiative  transport  equation  has  a  simple,  closed-form  solution 

/  =  — (  \-e~K'x)  (119) 

k'  v  ’ 

This  result  is  fundamental  to  the  well-known  radiative  transport  solution  method 
discussed  below. 

Tangent-Slab  Approximation 

The  method  of  solving  the  radiative  transport  equation  used  in  this  research  is  the 
well-known  tangent-slab  approximation  (Modest,  2003).  This  method  of  approximation 
splits  up  the  solution  into  to  equal  directional  components  (one  forward  and  on  reverse) 
along  a  particular  ray.  The  radiative  flux  is  calculated  by  assuming  that  each  grid  volume 
constitutes  a  thermodynamically  homogeneous  layer  (Greendyke,  1992),  whereby  the 
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spectral  radiative  flux  may  be  posed  in  terms  of  the  known  intensity  of  the  preceding 
layer 


(120) 


The  radiative  flux  in  a  given  layer  can  then  be  integrated  over  a  27i  steradians  solid  angle 
and  over  the  wavelength  interval  of  the  spectral  segment  under  consideration  to  obtain 
the  incident  radiative  flux.  This  integration  is  continued  in  marching  fashion,  from  the 
shock  to  the  wall,  along  a  path  normal  to  the  body. 

Conservation  Relation  for  Radiative  Energy. 

Finally,  the  radiative  transport  equation  may  be  more  generally  stated  in  a 
conservative  integral  form  as 

\\lv  (S  ■  n)dTdQ.  =  $  j(ev-K'vIv  )dVdCL 

n< r  n>v  (121) 

A  number  of  radiative  transport  schemes  which  are  more  spatially  and  directional  general 

than  the  two  so  far  discussed  may  be  developed  from  this  conservation  statement.  One 

such  approach  is  the  finite  volume  method  of  radiative  transport.  The  details  of  this 

method  and  its  application  in  the  present  investigation  are  briefly  discussed  in  the 

following  chapter  (Modest,  2003). 
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IV.  Methodology 


Flowfield  Solution  Procedure 

This  section  describes  the  implementation  of  the  foregoing  theory  in  a  suitable 
numerical  framework  via  the  flowfield  solver  NH7AIR.  The  system  of  equations  set 
forward  in  the  previous  chapter  (consisting  of  the  Navier-Stokes  equations  and  the 
various  source  terms,  therein  considered)  completely  describes  the  physics  of  a  flowfield 
in  thermochemical  nonequilibrium — within,  of  course,  the  inherent  limitations  of  the 
assumed  physical  models.  This  set  of  equations  however  is  not  suitable,  in  its  present 
form,  for  obtaining  numerical  solutions.  Therefore,  it  has  been  recast,  according  to  the 
method  presented  by  (Walters,  et  al.,  1990),  wherein  the  flowfield  equations  are 
numerically  solved  via  a  finite  volume  implementation  of  a  Roe-approximate  Riemann 
solver.  This  approach  involves  the  solution  of  the  local  Riemann  problem  at  the  cell 
interfaces  between  finite  volumes.  The  scheme  developed  for  perfect  gases  developed  by 
Roe  (1986)  has  been  extended  in  order  to  consider  thermodynamic  and  chemical 
nonequilibrium  in  three  dimensional  flows.  The  treatment  of  nonequilibrium  proposed 
by  Walters  and  his  colleagues  is  presented  here,  following  a  summary  of  the  method 
originally  proposed  by  Roe. 

Roe  Flux-Difference  Splitting. 

This  method  begins  by  casting  the  overall  flowfield  solution  in  terms  of  an 
ensemble  of  Riemann  problems  at  the  interfaces  between  the  finite  volume  cells  in  the 
solution  domain.  With  the  problem  thus  defined,  Roe  observes  that  -the  Riemann 
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solution  for  any  set  of  linear  conservation  laws  is  easily  computed.”  Therefore,  the 
derivation  of  a  suitable  scheme  begins  by  considering  the  linear  system 


<iw  .  <iw 

- +  A - =  0 

dt  dx 


(122) 


where  w  =  {p,pu,vp,...,pe}  is  the  vector  of  conserved  variables  and  A  is  the  constant 
Jacobian  matrix  defined  by  <3F  /  dw .  If  the  conserved  variables  to  the  left  and  right  of 
the  cell  interface,  wL  and  wR ,  are  known,  the  flux  difference  may  be  uniquely  expressed 
as 


Fr  -Fl  =  cxk  Ak t\  ,  (123) 

where  the  set  {ek}  contains  the  right  eigenvectors  of  A.  In  this  way,  the  contribution  of 
the  k-th  wave  to  the  flux  difference  in  given  in  terms  of  the  wave  strength  ak  and  wave 

speed  Xk .  It  is  evident  at  this  point  that  the  flux  at  the  cell  interface  (i+1/2)  may  be 
computed  by  either  expression 

If  1 2  (wz,  wfi)  =  Fl  akAkek  (124  a.) 


or 

F,:+i/2  (wi  > )  =  Fr  +  Z<+>  «/Ae,  •  (124  a.) 

By  averaging  the  two  foregoing  expressions 

®j'+i/2  (wi’ ) =  ^(Fr  +fl)  +  2  I \  I e<-  (124  a.) 

In  order  to  apply  the  foregoing  expression  to  a  nonlinear  problem,  one  must  first 
define  a  local  linearization  by  utilizing  A(wl,wr),  wherein  the  eigenvalues  and 
eigenvectors  of  the  resulting  linearization  not  only  satisfy  equation  (122)  but  also  the 
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eigenvectors  form  a  suitable  basis  set  whereby  the  Jump‘  in  the  conserved  variables 
across  the  cell  face  may  be  specified  by  the  linear  combination 

wr-wl=I«A  (125) 

k 

Ingeniously,  this  expression  returns  the  exact  solution  whenever  w£  and  wR  lie  on 

opposite  sides  of  a  flow  discontinuity.  Here  then  one  must  also  require  that  the  Rankin- 
Hugoniot  relationship  hold,  which  is 

Fr-Fl=S(wr-wl),  (126) 

where  S  is  the  shock  speed.  It  is  also  required  that  for  all  k 

Sak=  \ak  (127) 

This  statement  requires  that  all  ak  except  one  must  vanish.  Expressions  for  the  Roe- 
averaged  values  ak ,  Xk  and  ek  are  given  by  Roe  (1981).  Substitution  of  these  values 
back  into  the  final  expression  for  Fi+1/2 ,  gives  the  desired  solution  to  the  locally  linearized 

cell-interface  problem.  The  method  does  not  however  allow  for  the  finite  spatial 
distribution  of  expansion  wave  phenomena.  These  phenomena  can  be  accommodated  by 
an  entropy  fix  which  will  be  discussed  later  in  this  section. 

Having  briefly  reviewed  the  Roe  flux-splitting  scheme  developed  for  the  Euler 
equations,  consider  the  extension  of  this  scheme  to  accommodate  thermochemical 
nonequilibrium,  as  previously  discussed  (Walters,  et  al.,  1990).  The  governing  equation 
may  be  written  in  a  conservative  vector  from  in  2D  Cartesian  coordinates  as 

f <128> 


54 


where  Q  is  the  vector  of  conserved  variables,  W  is  the  vector  of  source  terms,  F  and  G 

/y  /v 

are  the  inviscid  flux  vectors  and  Fv  and  Gv  are  the  viscous  flux  vectors.  The  vectors  of 
conserved  variables  and  source  terms  are  given  below. 


'  A  ' 

'  d>x  ^ 

Pi 

d>2 

Pn 

d>N 

pu 

ZsN.z.E‘ 

pv 

,  W  = 

[ ZsN,z,e? 

pw 

Y.eKZ.E1 

P\en\ 

ti>\Dx  +  QT-y  i  +  Qe-v,  1  Qrad-V,  1 

..  ^ 

®eees  +  X  ^es,  (a),,.  +  Pele  +  Qr~e  ~  X  Qe-V 

s,r 

l  Pe o  J 

V  0  j 

(129) 


Before  introducing  the  inviscid  flux  vectors  it  is  necessary  to  introduce  some 
nomenclature.  The  arithmetic  average  of  a  quantity  /  is  calculated  from  the  left  and 
right  states,  as  indicated  by  subscripts  /  and  r ,  and  is  denoted  by  angled  brackets  below. 
Squared  brackets  denote  the  jump  of  quantity  /  across  the  cell  interface. 


</) 


(130  a.) 


f  =fr~fi 


(130  b.) 


The  approximate  Riemann  solution  requires  the  determination  of  the  cell  interface  fluxes. 
This  flux  is  calculated  as  a  summation  over  the  absolute  values  of  the  wave  speeds  A,  B 
and  C. 
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(131  a.) 


^(W) 


where 


+ 

LL  _U  LL  AlA  IL  _U5 


+ 


JiC 


(131  b.) 


The 


F 


vector  corresponds  to  the  eigenvalue  Aa=u  and  is  calculated  as 


indicated  below 


j 


r  P  A 


f  A  ] 

r  pjp  " 

Pi 

p2l  p 

Pn 

Pn1  P 

U 

u  u 

V 

+  J - -pu 

j 

V  ~Zy  U 

w 

w  -4  u 

en\ 

PPnJ  P 

enM 

pMeM  !  P 

Kk~a2 

0 

\  ^  J 

,  (132  a.) 


where 


M 

®  =  {u  U  +V  V  +W  W  U  +  'y1j 

- 1 

** 

- 1 

N 

-2>, 

— 1 

7=1 

p  J 

Z=1 

_p  l 

(132  b.) 


In  the  same  way  the  vectors  F  may  be  calculated  from  the  eigenvalues  AB  c  =  u±a 


B,C 
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=  M_L(  p  ± pa  u  )(«  +  «) 


1-  Me  J  2 a 


'  A 
A 


Pn 

u±£xa 

v±£ya 
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The  needed  Roe-averaged  quantities  were  calculated  as  indicated  below 
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Nonequilibrium  energy  terms  are  calculated  according  to 

[(pj%  I  p)4p 


en,  = 
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i  =  . 


(133) 


(134  a.) 

(134  b.) 

(134  c.) 


(134  d.) 


And  the  additional  thermodynamic  properties  of  enthalpy  and  entropy  are  calculated  as 
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(134  e.) 


where 
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(Tjp)  (e^) 


(134  f.) 


Finally  the  average  local  speed  of  sound  is  calculated  from  the  following  relation 


=  (/-!) 
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(134  g.) 


where 


(134  i.) 


and 


(134  j.) 

The  form  of  the  //  -direction  inviscid  flux  vector  G  may  be  found  after  the  same  manner 
as  the  £  -direction  flux  vector  developed  above.  Utilizing  a  thin  shear  layer 

approximation,  the  viscous  stress  tensor  Gv  can  be  written  as 
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where 
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(135  a.) 


(135  b.) 


g  =  w2  +  v2  +  w2 ,  (135  c.) 

and 

A  =  fjxu„  +V yv„  +  Vz^v •  (135  d.) 

MUSCL  Extrapolation. 

Second-order  spatial  accuracy  in  the  above  scheme  is  achieved  by  application  of 
the  MUSCL  extrapolation  (Van  Leer,  1979)  with  a  minmod  limiter  (Yee,  1987).  In 
essence,  the  MUSCL  extrapolation  replaces  the  piecewise  constant  interpolant  with  one 
which  is  piecewise  linear  on  the  solution  domain,  thus  increasing  the  solution  accuracy 
from  first-order  to  second-order.  However,  so  that  the  scheme  might  maintain  the 
property  of  being  total  variation  diminishing,  it  is  necessary  to  drop  to  first-order 
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accuracy  in  the  immediate  vicinity  of  a  flow  discontinuity.  This  limiting  procedure  is 
accomplished  by  application  of  the  minmod  flux  limiter  which  is  defined  by  the  function 

min  mod(x,  y)  =  sgn(x)  •  max  jo,  min  [|x|] ,  y  •  sgn(x)}  (136) 

One  may  thus  identify  the  slope  of  the  linear  interpolants  about  the  (i+1/2)  cell 
interface  by 


A  3  =  min  mod 
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(137  a.) 


(137  b.) 


And  finally  the  cell-center  values  of  the  conserved  quantities  may  be  extrapolated  the  left 
and  right  sides  of  the  cell  interface 


1  3 

U  .  1  o  ^  •  3 

J+2  ^  J+2 

(138  a.) 
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(138  b.) 

Entropy  Fix. 

In  order  to  eliminate  entropy-violating  phenomena  from  the  steady-state  solution, 
the  entropy  correction  i//(T)  is  applied  to  the  flux  scheme.  This  entropy  condition  is 

enforced  as  in  Josyula  and  Shang  (1993),  whereby  the  eigenvalues  are  cut  off  according 
to  the  relation 
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(139) 


The  isotropic  and  anisotropic  formulas  for  determining  8t  are  used  in  the  body-normal 
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(140) 


8n=8xJ~'  [|u  •  V£|  +  |u  •  V  77I  +  (a  /  2)  ( +  V  77)] 


and  body-tangential  directions 
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where  -  |u  •  Vk\  +  a  | Vk\  and  the  parameter  8X  is  assigned  values  of  0.5  and  0.01  in 

the  body-tangential  and  body-normal  directions,  respectively. 

Predictor-Corrector  Method. 

Time  integration  is  by  the  predictor-corrector  method  of  MacCormack  (1985).  It 
is  second-order  and  is  implemented  for  the  flux-splitting  method  in  these  steps: 

1)  Predictor  step 
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(141  a.) 

(141  b.) 


2)  Corrector  step 


At/;;1  =-At 
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(142  a.) 


U-'=±{U-J+U~‘+AU-'} 


(142  b.) 


Boundary  Conditions. 

The  Roe  flux-splitting  scheme  together  and  the  explicit  MacCormack  predictor- 
corrector  method  both  allow  for  the  use  of  explicit  boundary  conditions.  The  different 
boundary  and  initial  conditions  are  discussed  briefly  in  the  following  paragraphs.  Since 
the  nonequilibrium  flow  solver  is  based  on  the  finite  volume  method,  ghost  cells  are  used 
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to  implement  the  various  boundary  conditions  which  follow.  The  flow  solver  must  be 
supplied  with  a  grid  which  has  the  ghost  cells  explicitly  built  into  it.  Or,  said  another 
way,  the  code  treats  the  cells  along  the  edge  of  the  grid  as  ghost  cells. 

When  starting  the  flow  solver  it  is  necessary,  due  to  the  time-marching  nature  of 
the  solution  method,  to  provide  the  code  with  an  initial  condition.  The  flow  solver 
accepts  either  the  solution  from  a  previous  run  or  the  following  user-specified  data 
through  the  use  of  an  input  file  which  is  read  at  execution:  Twaii ,  Mref,  Lref  Tinf,  and  Pinf. 
At  the  first  time  step,  the  entire  domain  is  initialized  according  to  these  reference 
quantities. 

Given  the  ghost  cell  implementation  discussed  above,  the  wall  is  said  to  be 
located  at  the  cell  interface  between  the  ghost  and  first  interior  cells.  As  such,  the  no-slip 
boundary  condition  is  implemented  by  =cancelling  out‘  the  velocity  components  of  the 
adjacent  interior  cell  such  that  the  vector  average  at  the  cell  wall  is  identically  zero. 

ug=-ua  (143  a.) 

vg=-va  (143  b.) 

The  pressure  boundary  condition  at  the  wall  is  implemented  in  a  rather  straight 
forward  manner  by  assuming  a  zero  pressure  gradient  through  the  boundary  layer  to  the 
wall.  This  condition  provides  a  means  for  specifying  pressure  in  an  expedient  manner,  so 
long  as  the  boundary  layer  thickness  is  much  smaller  than  the  radius  of  curvature  of  the 
body  (White,  2006). 

Pg=Pa  (144) 
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The  isothermal  wall  boundary  condition  is  implemented  by  requiring  that  the  flow 
field  temperature  at  the  wall  (i.e.,  first  cell  interface)  be  equal  to  the  user-specified  wall 
temperature.  This  calculation  is  accomplished  by  a  simple  linear  extrapolation 

T  =2T  -T  (145) 

g  w  a  V  / 


The  ghost  cell  values  of  the  species-specific  vibrational  temperatures  and  the 
electron  temperature  are  set  according  to  an  assumed  quasi-adiabatic  condition.  This 
approach  is  reasonable  considering  that  these  modes  of  energy  transfer  are  much  less 
efficient  at  removing  thermal  energy  from  the  wall  than  collisions  with  heavy  particles. 
The  quasi-adiabatic  condition  is  enforced  by  setting  the  ghost  cell  value  equal  to  the 
adjacent  cell  value 


I  =T 

i,g  i,a 


(146) 


The  final  wall  boundary  condition  is  a  matter  of  specifying  the  nature  of  the 
chemistry  at  the  wall.  A  non-catalytic  boundary  condition  is  specified  at  the  wall,  which 
assumes  that  the  gradients  of  the  mass  fractions  are  zero  at  the  cell  interfaces.  This 
boundary  condition  yields  reasonable  results,  although  it  does  not  account  for  real  surface 
chemistry  effects  such  as  recombination. 

The  outflow  boundary  condition  is  somewhat  challenging  to  implement  due  to  the 
mixed  nature  of  the  solution  as  it  interacts  with  this  boundary.  Along  this  boundary,  the 
velocity  goes  from  zero  at  the  wall,  passing  through  the  subsonic  range  within  the 
boundary  layer,  to  supersonic  in  the  shock  layer.  In  the  subsonic  region,  the  solution  no 
longer  hyperbolic;  rather  it  is  parabolic.  That  is  to  say,  the  solution  on  the  boundary 
exhibits  a  certain  dependence  on  the  solution  within  the  domain  due  to  solution 
characteristics  capable  of  propagating  upstream.  Therefore,  it  is  desirable  to  split  the 
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boundary  condition  in  terms  of  the  local  Mach  number.  Where  the  local  Mach  number  is 
greater  than  one,  it  is  permissible  to  simply  project  the  value  in  the  solution  domain  into 
the  ghost  cell  by  some  suitable  projection.  Conversely,  where  the  local  Mach  number  is 
less  than  one,  a  boundary  condition  which  takes  into  account  the  dependence  would  be 
utilized.  One  such  boundary  condition  is  the  characteristic  outflow  boundary  condition 
of  Hirsch  (1987). 

According  to  this  approach  the  values  of  the  flow  quantities  at  the  outflow 
boundary  may  be  set  as  follows 

Pb  =  Pd  ~  Poao  (udn*  +  vdny )  ( 1 47) 

Pb  =  Pd -~r( Pd  ~ Pb)  (148) 

ao 

and  velocity  components  as 

ug=ua-2(udnx+vdny)nx  (149  a.) 

-  2  (Udnb  +Vdny)ny  ( 1 49  a0 

where  nx  and  n  are  the  components  of  the  outward-facing  normal  vector  of  the  boundary 

cell.  The  NH7AIR  code  currently  implements  the  supersonic  outflow  boundary 
condition.  It  was  proposed  that  the  outflow  boundary  condition  be  improved  according  to 
the  subsonic  implementation  discussed  above.  This  implementation  was  briefly  pursued. 
However,  this  introduced  unexpected  numerical  instabilities  which  were  not  able  to  be 
resolved  in  a  timely  manner.  While  the  mixed  boundary  condition  would  have 
represented  the  physical  situation  more  accurately,  the  original  boundary  condition  was 
accepted  for  the  sake  of  moving  forward  with  more  central  research  tasks. 
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The  implementation  of  the  inflow  boundary  condition  is  relatively 
straightforward.  Since  the  governing  equation  set  is  hyperbolic  in  nature  ahead  of  the 
shock  layer,  the  inflow  boundary  condition  has  no  dependence  on  the  solution  within  the 
domain.  As  such,  it  is  only  required  that  the  values  along  this  boundary  (i.e.,  ghost  cell 
values)  be  specified  as  the  freestream  values. 


Grid  Adaptation. 

There  are  two  regions  of  great  interest  in  the  flow  surrounding  a  reentry  vehicle. 
First,  the  region  near  the  bow  shock  is  of  great  importance  due  to  the  relaxation  processes 
which  occur  just  downstream  of  it.  These  relaxation  processes  greatly  influence  the 
radiative  heating  predicted  by  the  radiative  transport  method.  The  second  region  of 
interest  is  the  boundary  layer,  which  is  known,  of  course,  to  determine  the  convective 
heating  predicted  by  the  flow  solver.  Therefore,  it  is  necessary  to  adapt  the  grid  to 
adequately  compute  these  and  other  important  quantities.  Grid  adaptation  is 
accomplished  upon  initiation  of  a  given  calculation  session,  according  to  the  method  of 
Gnoffo,  et  al.,  (1993)  and  as  implemented  in  the  NH7AIR  code  by  Komives  and 
Greendyke  (2009;  Komives  2009). 

The  algorithm  utilizes  four  user-specified  parameters  in  order  to  perform  the  grid 
adaptation  on  the  k  <  K  cell  faces  in  the  wall-normal  direction.  The  first  such  quantity  is 
the  cell  Reynolds  number,  N  Rece// ,  which  determines  the  first  cell  size  according  to 


A«(l)  = 


(150  a.) 


where 
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(150  b.) 


^Rece;/  = 


paAn 

F 


A  few  notes  are  needed  here  before  proceeding.  First  let  n{]}  ( k )  denote  the  original 
distance  between  the  body  surface  and  the  k-th  cell  center  along  the  coordinate  C,.  The 
value  h(k)  =  (k)ln®  (K)  denotes  the  nondimensional  distance  in  the  C,  direction. 

Furthermore,  let  Ah(k)  =  n (£  +  1/2)  —  k  — 1/2) define  the  width  of  the  k-th  cell  with 

the  interpretation  that  the  indices  k+1/2  and  k- 1/2  refer  to  the  outer  and  inner  cell  edges, 
respectively.  The  second  parameter,  F  bl ,  specifies  the  fraction  of  the  total  number  of 

cells  which  are  placed  in  the  boundary  layer  by  the  mapping,  Kbl  =  FblK .  The  following 


function  controls  the  cell  growth  in  the  mapping  of  cells  into  the  boundary  layer 


Ah[k) = 


1  +  C  sin 


(*-!)■ 


Ah[k-X) 


(151) 


where 


C  = 


An(l) 


-1 


(152) 


With  the  cell  individual  cell  widths  calculated  according  to  the  above  expression,  the 
distribution  of  h  is  obtained  via 


n(k  +  M2)  =  ^An(/)  (153) 

/= i 

This  transformation  provides  for  gradual  cell  growth  in  the  boundary  layer.  This  growth 
slows  down  as  the  edge  of  the  boundary  layer  is  approached,  such  that  the  remaining  cells 
past  k  >  Kbl  are  equally  spaced. 
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The  second  transformation  is  designed  to  resolve  the  bow  shock  by  grouping  cells 
near  the  h  =  Fsh  location.  First,  a  renormalization  is  preformed  to  fix  range  of  zeta  exactly 

between  zero  and  one.  This  normalization  is  done  by  dividing  each  cell  face  location  by 
the  cell  face  on  the  outer  boundary, 

h(k  + 1/2)/  n(K  + 1/  2) - >h(k  +  \l  2).  (154) 

With  the  renormalization  complete,  the  transformation  is  performed  according  to 

h(k  +  \l =  (\—e(k  +  \/2)h(k  +  \l 2^  +  Fshe(k  +  H 2)  (155  a.) 

where 

e{k  +  \l2)  =  h2{k  +  \l2)(\-h{k  +  \l2))e0  (155  b.) 

This  expression  introduces  the  fourth  user-specified  parameter,  ,  which  controls  how 

tightly  the  cells  are  group  about  the  shock.  This  parameter  must  be  chosen  with  care  to 
ensure  that  the  grid  does  note  fold  back  onto  itself  near  the  shock. 

The  final  transformation  returns  the  distribution  of  cell  centers  in  the  original 
dimensions  of  the  grid.  A  scaling  factor  is  used  to  ensure  that  the  captured  shock  lies  at 
the  specified  fraction  Fsh  of  the  distance  between  the  body  and  the  outer  boundary.  This 
transformation  is  performed  according  to 

(2)/  \  n^(*)n(k) 

n{2\k)  = - w  W  (156) 

FSh 

where  h1'1  (*)  is  the  location  on  the  original  grid  where  the  captured  shock  is  first  sensed. 
Finally,  interpolation  and  extrapolation  are  used  to  map  all  the  old  grid  points 
x[n^j  (k  + 1  / 2) j into  the  new  grid  x {n^j  (k  + 1  / 2) j . 
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Radiation  Transport  Solution  Procedure 


Two  methods  of  solving  the  radiative  transport  equation  were  developed  and  used 
to  obtain  results  in  this  research.  The  tangent  slab  method,  discussed  in  the  previous 
chapter,  was  used  to  obtain  radiative  solutions  both  coupled  and  uncoupled  with  the 
nonequilibrium  flow  solver.  The  results  obtained  by  the  finite  volume  method  for 
radiative  transport  were  so  obtained  in  an  uncoupled  fashion  but  utilizing  the  coupled 
flow  fields  resulting  from  coupling  the  two-flux  method  with  the  nonequilibrium  flow 
solver. 

Tangent  Slab  Method. 

The  implementation  of  the  tangent  slab  method  within  the  context  of  a  reentry 
shock  layer  is  rather  straightforward  and  is  a  stardard  part  of  spectroscopic  codes  like 
SPRADIAN.  The  expression  presented  in  equation  (119)  may  be  easily  evaluated  in  a 
marching  fashion  toward  and  away  from  the  body.  With  the  radiative  intensities  thus 
obtained,  the  radiative  source  terms  may  be  evaluated  and  coupled  with  the 
nonequilibrium  solver  as  discussed  in  the  following  section. 

Finite  Volume  Method. 

The  finite  volume  method  of  radiative  transport  (FVMR)  is  based  on  the 
conservation  relation  for  radiative  energy  given  in  Chapter  III,  and  its  development  for 
use  in  the  axisymmetric  flowfield,  as  part  the  present  research  effort,  is  a  unique  feature 
of  this  work. 

\\iv  {s  ■  h)dTdn  =\\(ev-K'viv  )dVdn 

o,r  n,y  (157) 
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Assuming  that  intensity  is  constant  across  a  given  face  of  the  element,  as  well  as  over  the 
solid  angle  D.i ,  in  the  i-th  direction,  the  foregoing  conservation  equation  can  be  restated 
in  the  discretized  form 

2 1,M  (5  ■ », )  A  =  (*„  - K  'Iy„ )yn,  (158  a.) 

where 

y  =  J  sdQ.  (158  b.) 

n, 

In  the  above  expressions,  the  unit  normal  vector  and  area  associated  with  the  k-th  cell 
face  are  given  by  nk  and  Ak ,  and  the  unit  normal  vector  and  total  solid  angle  vector 

associated  with  the  i-th  direction  are  given  by  s  and  si .  Figure  1  illustrates  the  method 

by  which  the  foregoing  unit  vectors  are  assigned  to  the  spatial  and  directional 
discretization  schemes.  The  intensities  at  the  face  centers  //„■  are  related  to  those  at  the 
volume  centers  according  to  the  step  scheme  where  a  positive  or  negative  dot  product 
(y.  •  hk )  indicates  a  flux  out  of  or  into  a  cell,  respectively.  First,  it  is  assumed  that  for 

intensities  leaving  the  control  volume  P  the  intensity  at  the  k-th  face  is  equal  that  of  the 
subject  volume ‘s  cell-center  intensity  Ipi  in  the  i  direction.  Then  for  intensities  entering 
the  subject  volume  it  is  possible  to  take  Iu  to  equal  the  Ipi  of  the  appropriate  neighboring 
cell. 
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Figure  1.  Spatial  and  Directional  Discretization  for  Finite  Volume  Method  for  Radiation 

With  the  spatial  and  directional  domains  specified  as  above,  it  is  possible  to  make 
the  appropriate  substitutions,  and  thus  solve  for  the  cell-center  intensities  7P;  explicitly  in 
terms  of  the  neighboring-cell  intensities  /* j. 
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Finally,  the  incidence  radiation  and  radiative  flux  are  calculated  according  to 


(159) 


and 


G  =  1,1  Q. 

vp  .  Vip  l 


Vvp^hipSi 


(160  a.) 


(160  b.) 


Suitable  boundary  conditions  for  the  radiative  transport  equation  can  be  easily  obtained  in 
a  similar  manner  and  are  stated  here,  where  -s-”  a  surface  quantity. 

(161) 
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Numerical  Solver  Details. 


Given  the  form  of  equation  (159),  it  is  possible  to  solve  for  the  cell-center 
intensities  /  .  utilizing  one  of  two  methods.  The  first  method  would  be  to  obtain  a 

solution  iteratively,  by  guessing  an  initial  solution  for  center  intensities  /  . .  Using  this 
initial  guess,  the  solution  is  marched  forward  using  the  cell-centered  intensities  /  .  at 
each  (n  - 1  )th  step  to  approximate  the  cell-neighbor  intensities  Ivkl  at  the  nth  step  and 

thus  obtain  the  updated  cell-center  intensities.  This  process  continues  until  the  solution 
reaches  an  acceptable  level  of  convergence. 

A  second  method  would  be  to  solve  the  linear  system  of  equations  created  by 
expressing  equation  (159)  in  matrix  form.  Where  the  cell-center  intensities  /  .  are  now 

considered  as  Iv(I  Jy ,  the  intensity  at  grid  point  (/,./) ,  and  the  neighbor-cell  intensities  Ivki 

are  taken  at  the  neighboring  grid  points  (/  + 1, J),  (/  - 1,  j\  (/,  j  + 1),  and  (/,  j  - 1) .  This  scheme 
results  in  the  banded  linear  system  with  diagonal  sub-matrices  defined  as  follows 
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Fl  +k'V  Q  Fl 

1 1 ,j  '  ^  v  i,F*i  1  i+i,j 
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pi  _j(nk*si),if  <  0 

u*u  i  Q 


(162  b.) 


F'u  =  Z 

k,if>  0 


(162  c.) 
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This  system  has  a  bandwidth  which  corresponds  to  the  number  of  J  grid  points  in 
each  grid  row  I.  The  elements  F‘k  correspond  to  the  inner  product  of  the  area  normal 

vector  hk  with  the  transmission  direction  vector  s, ,  where  F:  7  is  calculated  based  on  the 

outgoing  radiation  and  all  off-diagonal  Flk  elements  are  based  on  the  incoming  radiation. 

The  solver  used  in  this  code  utilized  Gaussian  elimination  with  pivoting  to  improve 
numerical  stability. 

Additional  Geometrical  Considerations. 

One  of  the  additional  considerations  for  implementing  the  present  transport 
scheme  is  the  method  by  which  the  cell-center  field  of  view  is  discretized  by  specification 
of  both  the  magnitude  and  direction  of  the  differential  solid  angle  fX .  For  the  purposes 

of  this  dissertation,  the  following  assumptions  have  been  made  regarding  £X . 

1.  Each  differential  solid  angle  fX  is  defined  by  a  vector  Si  which  specifies  its 
orientation  and  angular  extent. 

2.  Each  Si  occurs  at  the  centroid  of  the  differential  solid  angle 

3.  Each  fX  consists  of  a  continuous,  regular”  solid  angle 

4.  The  topology  of  Q;  upon  the  unit  sphere  is  such  the  distance  between  the 
edge  dQ.i  and  the  centroid  of  Q,  is  minimized. 

5.  Each  Q;  contributes  Ak  /  N  steradians  to  the  total  FOV,  where  N  is  the 
number  of  differential  solid  angles. 

6.  The  FOV  should  be  such  that  the  resulting  Si ‘s  are  spaced  at  regular  intervals 
within  the  FOV 
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Generally  speaking,  the  finite  volume  method  radiative  transport  scheme  poses  no 
restriction  on  the  number  of  transmission  direction  which  may  be  considered.  For  the 
purposes  of  this  dissertation,  six  transmissions  have  been  considered — one  each  in  the 
positive  and  negative  z,  r  and  9  directions.  This  is  the  minimum  number  needed  to 
investigate  possible  3D  effects.  This  selection  of  transmission  directions  is  also 
advantageous  for  comparisons  with  experimental  data  which  are  most  often  collected  at 
orientations  which  are  normal  to  coordinate  system  basis  vectors. 

One  particularly  challenging  aspect  of  implementing  the  FVMR  scheme  in  an 
axisymmetric  coordinate  system  is  related  to  propagating  the  solution  across  the  grid 
singularity  created  by  line  of  symmetry.  The  problem  lies  in  two  areas.  First,  there  is  a 
geometric  constriction/relief  effect  in  the  r-direction  and  second,  the  simulation  of 
adjacent  neighbor  cells  in  the  9  direction.  Because  of  symmetry,  it  is  desirable  to 
compute  the  solution  on  a  wedge-shaped  region  of  the  total  flow  field.  Therefore,  the  in¬ 
direction  faces  become  increasingly  small  and  vanish  as  the  centerline  is  approached. 
Conversely,  the  r-direction  faces  become  increasingly  large  moving  away  from  the 
centerline.  Here  the  difficulty  becomes  apparent  as  two  geometric  effects  begin  to 
influence  the  solution.  Fluxes  which  approach  the  singularity  are  constricted.  Because 
the  FVMR  conserves  radiative  energy,  the  intensity  exiting  through  the  shrinking  r-faces 
will  become  increasingly  large  and  even  become  unbounded  at  the  singularity.  As  fluxes 
depart  from  the  centerline,  there  is  a  geometric  dissipation  effect. 

The  contribution  of  radiative  flux  from  these  adjacent  cells  is  approximated  by 
assuming  that  the  local  intensity  at  the  9  -direction  cell  face  is  equal  to  the  intensity  at  IJ 
cell  center  and  acting  through  the  projection  of  the  9  faces  onto  the  r-z  plane.  However, 
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a  straightforward  application  of  the  step  scheme  above  would  not  work  here  since  the 
neighboring  cells  in  the  9  direction  are  by  symmetry  the  same  as  the  cell  under 
consideration.  Considering  the  negative  r  transmission  direction  first,  it  is  noted  that  the 
dot  product  of  the  negative  r-direction  unit  vector  with  the  unit  vector  of  the  9  faces 
(sr_  -  fig)  is  negative,  and,  therefore,  the  associated  Fr ,  matrix  element  in  equation  (162) 

preserves  the  sign  convention  already  establish  by  the  step  scheme.  However,  the 
situation  is  more  difficult  in  the  positive-r  transmission  direction  first.  In  this  case,  the 
dot  product  (sr+  •«„)<(),  and  so  there  is  an  influx  of  radiative  flux  from  the  IJ  cell  center. 

Unfortunately,  it  would  violate  the  sign  convention  of  the  step  scheme  to  assign  a 
negative  value  of  FI ,  to  the  IJ  diagonal.  Thus,  the  following  alternative  approach  is 

proposed.  The  approach  begins  with  the  observation  that  the  on-diagonal  elements  of 
equation  (162)  consist  of  two  terms:  the  cell-face  view  factor  FtJ  and  the  cell-center 

value  of  absorption  k'Vi  yQ .  Furthermore,  it  is  required  that  under  the  step  scheme 
Fjj>  0  and  generally  k'V/  /Q.  >  0 .  While  it  is  not  permissible  to  assign  the  negative 
value  (y.+  •  he )  <  0  to  the  on  diagonal  view  factor  Ff , ,  it  is  permissible  to  interpret  this 

incoming  radiative  flux  as  a  being  a  contribution  to  the  absorption  expressed  by  the  term 
/c'VjjQ..  Thus,  it  proposed  to  model  this  absorption  via  a  geometric  absorption  term 

which  is  equal  in  magnitude  and  opposite  in  sign  to  (sr+  ■  h0  ) . 

Further  difficulty  arises  in  the  case  of  radiative  fluxes  in  the  9  direction  due  to  the 
fact  that  only  a  one  cell  wide,  wedge-shaped  region  of  the  flow  field  is  modeled.  Given 
the  above  considerations,  future  implementation  of  the  FVMR  scheme  should  avoid  an 
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axisymmetric  grid  implementation  such  as  the  one  described,  in  favor  of  a  full  3D  grid 
with  no  singularity. 

Radiative-Gasdynamic  Coupling 

In  this  section  a  suitable  coupling  scheme  between  the  radiation  and  gasdynamic 
solutions  is  developed.  Pai  (1966)  offers  a  complete  treatment  of  fully-coupled  radiation- 
gasdynamic  system;  however,  a  loosely-coupled  scheme  is  used.  In  such  a  scheme, 
radiation-related  source  terms  are  static  (or  ramped)  over  several  iterations  of  the  flow 
solver  and  updated  according  to  user-specified  criteria,  such  as  order  of  convergence  or 
number  of  iterations.  This  approach  is  allowable  since  time-accurate  simulations  are  not 
being  pursued.  However,  under  this  approach,  convergence  is  not  guaranteed  and 
solutions  may  not  necessarily  be  unique.  Given  the  very  stiff  nature  of  the  radiation 
terms,  numerical  challenges  are  likely  to  be  encountered  in  such  a  coupled  situation, 
especially  when  using  shock-capturing  techniques  (Gnoffo,  et  al.,  2009). 

The  additional  consideration  of  radiation  in  a  flowfield  requires  the  tracking  of  a 
new  energy  transfer  mechanism,  namely  that  energy  which  is  transported  through  the 
solution  domain  due  to  radiation.  In  order  to  have  a  coupled  scheme,  that  energy  must 
show  up  in  the  flowfield  equations.  In  practice,  one  is  interested  in  the  total  radiative 
source  term  in  addition  to  the  radiation-electronic  energy  source  term  and  radiation- 
species  vibrational  energy  source  terms,  since  these  are  the  ones  which  are  needed  to 
couple  the  radiation  transport  and  nonequilibrium  flow  solver  codes. 

Recalling  the  conservation  relation  for  radiative  energy,  the  total  amount  energy 
lost  from  the  flow  field  due  to  radiation  may  be  computed  by  integrating  the  divergence 
of  the  spectral  intensity 
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(163) 


n,  r  v 

Alternately,  the  net  emission  may  be  integrated 

oo 

a^^llS^-K'I^dVdSi  (164) 

nv  o 

This  energy  represents  the  difference  between  emitted  and  absorbed  radiation.  Thus,  if 
the  fluid  emits  more  radiation  than  it  absorbs,  then  the  energy  is  lost  in  a  phenomena 
known  as  radiation  cooling”.  Conversely,  the  fluid  may  be  heated  by  absorbing 
radiation.  In  many  situations,  the  absorption  of  radiation  is  considered  to  be  negligible, 
relative  to  other  effects.  When  this  is  the  case,  the  media  is  said  to  be  optically  thin,  and 
it  is  possible  to  evaluate  the  source  term  solely  in  terms  of  the  emission  coefficient  via  the 
simple  algebraic  formula 

Qtotal=A™vV  (165) 

Coupling  with  the  vibrational  energy  equation  occurs  through  the  radiative 
species-vibrational  source  terms 

oo 

Qvib,s  =  J  J  J  «...  (£v,s  -  Kv,s  \s  )dvdVd£l  (166) 

nr  o 

where  svs  and  k  s  denote  the  contribution  of  species  s  to  the  emission  and  absorption 
coefficients  and  av  s  represents  proportion  of  energy  which  is  responsible  for  exciting  or 

damping.  This  calculation  may  be  accomplished  by  considering  the  separability  of 
internal  energy,  whereby  a  portion  of  the  energy  hv  contained  in  the  emitted/absorbed 

photons  can  be  attributed  to  the  change  in  vibrational  energy  resulting  from  a  given 
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transition.  Consider  that  the  energy  at  a  particular  energy  level  which  is  the  sum  of  the 
electronic,  vibrational  and  rotational  terms  as  shown 

E(e,v,J)  =  T'+G,+F,.  (167) 

Thus  the  energy  which  must  emitted  or  absorbed  by  a  photon  in  order  to  undergo  a  state 
transition  is  given  as 

hv  =  AE  =  (T, -T,)  +  (G,  -  G,)+(Fj. -F,)  (168) 


It  then  follows  that  the  proportion  of  that  energy  which  contributes  to  the  excitation  of  the 
vibrational  mode  is  given  as 


avib  = 


{g,-gv) 


(Te-Te)  +  (Gv-Gv)  +  (Fr-Fj) 


(169) 
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V.  Computer  Implementation 


Having  discussed  the  necessary  background  theory  and  various  methodological 
aspects  of  this  research  effort  in  the  previous  two  chapters,  the  present  chapter  outlines 
how  these  concepts  have  been  implemented  in  a  computer  code.  In  the  course  of 
outlining  this  implementation  process,  the  basic  structure  and  function  of  the  two  baseline 
codes  are  discussed,  and  the  manner  in  which  they  have  been  modified  and  coupled  is 
explained.  The  basic  flow  of  the  resulting  computer  program  is  illustrated  below  in 
Figure  2. 


Figure  2.  Flowchart  for  Overall  Computer  Program 

This  chapter  is  organized  into  three  sections:  two  sections  correspond  to  the 
modifications  made  to  the  baseline  flow  field  and  spectroscopic  solvers,  NH7AIR  and 
SPRADIAN,  respectively,  and  another  section  which  follows  the  development  of  a 
radiation  solver  which  calculates  the  solution  to  the  radiative  transport  equation  and 
handles  the  passage  of  data  between  the  flow  solver  and  spectroscopic  code.  The  various 
elements  of  the  computer  program  are  listed  below  in  Table  2. 
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Table  2.  Summary  of  Computer  Program  Elements;  highlighting  indicates  those  aspects 
of  the  computer  code  which  have  been  modified. 


Flow  Solver  (NH7AIR) 

Radiation  Solver 

Spectroscopic  Solver  (SPRADIAN) 

subroutine  a360 

module  rad_parameters 

module  size_def_mod 

subroutine  alignshock 

module  rad_vars 

module  struct_def_mod 

subroutine  be 

subroutine  rad_solver 

module  interface_mod 

subroutine  const 

subroutine  rad_TS 

subroutine  radipac 

subroutine  cvmgp 

subroutine  rad_FVM 

subroutine  emis_absb 

subroutine  datin 

subroutine  rad_couple 

subroutine  atom_bb 

subroutine  estdt 

subroutine  simpson 

subroutine  atom_bf 

subroutine  fsiroe 

rad_comon 

subroutine  atom_ff 

subroutine  fsjroe 

rad_parameters 

subroutine  atom_noneq 

subroutine  gcomon 

rad_solver90 

subroutine  diatom_bb 

subroutine  gmtry 

radipac6X90 

subroutine  calc_diatom_dist 

subroutine  gridin 

subroutine  calc_diatomic_bb 

subroutine  hisstr 

subroutine  diat_eimp_exct 

subroutine  iviserg 

subroutine  cros_ab 

subroutine  jviserg 

subroutine  diatom_noneq 

subroutine  1 

subroutine  diatom_read 

subroutine  Imitri 

subroutine  H_bb 

subroutine  Imitrj 

subroutine  intpll 

program  main 

subroutine  minv 

subroutine  p3dwr 

subroutine  monatom_read 

subroutine  parse 

subroutine  simp 

subroutine  plotc 

subroutine  taint 

subroutine  sourcet 

subroutine  tri_cont 

subroutine  stvar 

subroutine  triatom_read 

subroutine  transp 

subroutine  vuv_bf 

subroutine  vtkio 

subroutine  wrstte 

subroutine  wrtout 

Flow  Solver 

The  discussion  of  the  major  segments  of  the  developed  computer  code  begins 
with  the  nonequilibrium  flow  solver  NH7AIR,  since  this  portion  of  the  code  most  directly 
controls  the  overall  execution  of  the  program.  From  the  perspective  of  the  solution 
methodology,  this  aspect  of  the  resulting  code  is  not  surprising.  Consider  the  nature  of 
the  coupling  between  the  flow  field  and  the  radiative  solutions.  At  the  length  and  time 
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scales  involved  when  considering  reentry  situations,  the  radiation  solution  exhibits  an 
elliptic  behavior  in  the  solution  domain.  Furthermore,  the  radiation  solution  may  be 
assumed  to  update  =instantly‘  to  new  flow  field  conditions;  whereas  several  iterations  of 
the  flow  solver  are  required  to  allow  the  flow  field  to  adjust  to  new  radiative  source 
terms.  The  computer  code  which  results  from  this  methodology  is  one  wherein  the 
radiation  solver  and  spectroscopic  solver  function  as  subroutines  which  are  periodically 
called  by  the  main  program — the  flow  solver — in  order  to  update  the  radiative  source 
terms.  A  brief  description  of  the  most  significant  aspects  of  the  baseline  flow  solver 
follows  in  order  to  facilitate  the  subsequent  discussion  of  modifications  made  thereto. 
Figure  3  contains  a  flow  chart  which  illustrates  the  logical  arrangement  of  the  most 
important  subroutines  in  the  flow  solver. 

The  main  program  is  contained  in  the  Fortran  file  so  named  main./.  As  is 
customary,  the  main  program  coordinates  the  execution  of  the  overall  computer  code  by 
performing  the  primary  input/output  functions  and  calling  the  various  subroutines 
contained  in  the  program.  The  subroutine  datin  is  the  first  called  and  is  responsible  for 
reading  the  input  deck  and  the  restart  files,  consisting  of  a  grid  and  solution  files  from 
previous  runs.  Initial  and  boundary  conditions  are  supplied  by  the  subroutine  be,  which 
is  called  at  restart  and  at  each  time  step.  The  main  program  loop  consists  of  calls  to  the 
subroutines  indicated  in  Figure  3. 


80 


Figure  3.  Flowchart  for  Nonequilibrium  Flow  Solver  (NH7AIR);  bold  indicates  areas  of 

the  code  affected  by  modifications. 


Most  of  the  subroutines  listed  in  the  main  loop  have  names  which  suggest  the  function(s) 
they  perform.  The  subroutine  estdt  estimates  the  local  time  step  based  on  the  CFL 
criteria  specified  by  the  input  file.  Transport  coefficients  for  heat  transfer  and  viscosity 
are  calculated  by  transp.  The  subroutine  L  has  a  rather  opaque  name  yet  is  critically 
important.  This  subroutine  accomplishes  the  time-integration  of  the  solution  and 
performs  calls  to  several  other  supporting  subroutines.  The  calculation  of  convective 
fluxes,  via  the  Roe-averaging  method  discussed  in  the  previous  chapter,  is  performed  by 
fsiroe  and  fsjroe,  where  the  subroutines  Imtri  and  limtrj  implement  a  minmod  limiter  as 
described  previously.  Calculation  of  the  viscous  fluxes  is  performed  by  iviscrg  and 
jviscrg.  Subroutine  stvar  backs  out  the  state  variables  from  the  vector  of  conserved 
quantities  and  sourcet  calculates  the  various  source  terms  used  in  L  to  update  the  flow 
field  solution.  Finally,  various  subroutines  write  restart  files  and  desired  output  files. 
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Modifications. 

The  modifications  made  to  the  baseline  flow  solver  were  primarily  concerned 
with  input/output  functions,  the  passage  of  data  between  the  flow  and  radiation  solvers, 
and  the  introduction  of  radiative  source  terms  into  the  flow  solver.  Listings  of  these 
modifications  may  be  found  in  Appendix  A  and  are  summarized  below  in  Table  3. 


Table  3.  Summary  of  Modifications  Made  to  Baseline  Flow  Solver. 


Description 

Subroutine 

Appendix  Entry 

Changes  to  input  file 

datain 

Table  11 

Radiation  restart  (Read) 

datain 

Call  radiation  solver 

main 

Table  12 

New  output  files 

main 

Modifications  to  vtk  output 

vtkio 

Radiation  restart  (Write) 

wrtout 

Here  follows  a  brief  discussion  of  these  modifications.  The  first  set  of  modifications 
affect  the  subroutine  datain,  wherein  additional  read  statements  were  needed  in  order  to 
input  solution  parameters  associated  with  the  setup  and  execution  of  the  radiation  solver. 
This  subroutine  was  further  modified  in  order  to  read  radiation  restart  data  into  the 
program.  Radiation  restart  files  are  written  by  statements  added  to  the  subroutine  wrtout. 
In  addition,  various  modifications  were  made  to  main  and  the  subroutine  vtkio  in  order  to 
output  quantities  of  interest  associated  with  the  radiation  solution.  The  most  significant 
modifications  were  to  add  a  call  to  the  radiation  solver  in  the  main  loop  and  to  update  the 
source  terms  in  order  to  account  for  the  effects  of  radiation.  These  modifications  were 
made  to  the  main  program  and  to  the  subroutine  sourcet. 
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Radiation  Solver 


The  discussion  of  the  implementation  effort  now  turns  to  the  development  of  the 
radiation  solver,  which  is  responsible  for  solving  the  radiative  transport  equation  (RTE) 
and  updating  the  radiative  source  terms  in  the  energy  conservation  equations  of  the  flow 
solver.  The  radiation  solver  is  called  from  within  the  main  loop  of  the  flow  solver  and 
thus  is  subordinate  in  program  hierarchy  to  the  flow  solver.  The  sequence  of  the 
subroutine  components  (and  their  interaction  with  program  segments  outside  of  the 
radiation  solver)  is  illustrated  below  in  Figure  4. 


Figure  4.  Flowchart  for  Radiation  Solver  ( radjsolver ) 


As  illustrated  in  the  preceding  flowchart,  the  execution  of  the  radiation  solver 
subroutine  consists  in  the  sequential  performance  of  the  following  functions:  setup,  solve 
RTE,  and  calculate  source  terms.  This  functional  delineation  provides  a  convenient 
framework  for  discussing  the  subcomponents  of  the  radiation  solver. 
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The  first  function  to  be  performed  is  to  setup  the  radiation  solver  in  terms  of 
various  user-specified  parameters  and  flow  field  variables — all  of  which  are  passed  by 
the  main  program  as  arguments  to  the  subroutine  rad_solver.  Receiving  these  values 
from  the  main  program,  the  remainder  of  the  setup  function  is  carried  out  by  the  modules 
rad_vars  and  rad _parameters  and  by  the  subroutine  specify _parameters.  The  module 
rad_vars  declares  a  variety  of  variables  which  are  shared  by  the  local  subroutines 
contained  in  rad_solver  and  used  in  the  calculation  of  the  radiative  intensities  and  source 
terms.  Similarly,  the  module  rad parameter  declares  various  parameters  associated  the 
solution  procedure;  most  importantly,  it  assigns  character  strings  to  the  arrays  which 
identify  the  various  radiative  mechanisms  used  in  both  the  solution  of  the  radiative 
transport  solution  procedure  and  by  the  spectroscopic  solver.  The  subroutine 
specify parameters  determines  the  number  of  active  radiative  mechanisms  and  parses  the 
above  string  arrays  in  order  to  make  them  more  amenable  to  the  format  of  the  data 
structures  in  the  subroutine  rad_solver. 

Subsequent  to  the  setup  of  the  radiation  solver,  the  subroutine  rad_solver  calls 
one  of  two  subroutines  in  order  to  solve  the  radiative  transport  equations:  rad_TS  or 
rad_FVM.  These  two  subroutines  are  addressed  now  in  turn.  The  first  subroutine, 
rad_TS,  calculates  the  body-normal  components  of  the  radiative  intensity  in  the  forward 
and  reverse  directions.  The  general  flow  of  this  subroutine  is  summarized  below  in 
narrative  and  also  in  the  pseudo-code  provided  Tables  17-22  of  Appendix  A.  Upon 
declaring  and  initializing  variables  which  are  local  to  rad_TS,  a  set  of  nested  loops  is 
executed  in  order  to  calculate  radiative  intensity.  The  outermost  loop  advances  the 
solution  procedure  from  one  body-normal  path  to  the  next.  At  the  next  level  down,  the 
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second  loop  selects  the  direction  of  integration — namely,  the  forward  and  reverse 
directions  along  the  path  specified  by  the  outer  loop.  The  inner  loop  advances  the 
solution  from  one  point  to  another  along  the  integration  path  and  in  the  direction 
specified  in  the  outer  two  loops.  At  each  point  along  the  integration  path,  it  is  necessary 
to  obtain  the  values  of  the  spectral  emission  and  absorption  coefficients  from  the 
spectroscopic  subroutine  radipac,  and  then  to  obtain  the  radiative  intensity  from  the 
tangent  slab  solution  to  the  radiative  transport  equation.  Note  that  the  calculation  of  the 
emission  and  absorption  coefficients  is  relatively  expensive  and  is  only  performed  only 
once  at  each  point  in  the  flow  domain  for  a  given  iteration  of  the  radiation  solver. 
Finally,  the  radiative  source  terms  are  calculated  from  the  local  emission,  taking  into 
account  the  absorption  of  radiation  from  the  two  transmission  directions  considered. 

The  second  subroutine,  rad_FVM,  calculates  the  components  of  the  radiative 
intensity  in  the  positive  and  negative  coordinate  directions  of  a  cylindrical  coordinate 
system  aligned  with  the  centerline  of  the  flow  field.  These  directions  were  chosen  for 
convenience.  The  general  flow  of  rad_FVM  is  summarized  below  in  narrative  and  also 
in  the  pseudo-code  provided  Tables  23-27  of  Appendix  A.  As  before,  the  subroutine 
begins  by  declaring  and  initializing  the  local  variables  needed  by  rad_FVM.  The  next 
step  is  to  calculate  the  emission  and  absorption  coefficients.  The  procedure  for  solving 
the  radiative  transport  equation  is  fundamentally  different  for  the  FVMR  scheme  and  the 
tangent  slab  method.  Whereas  the  tangent  slab  solution  procedure  utilizes  an  analytical 
approximation  in  a  local  integration,  the  FVMR  solves  the  radiative  transport  equation 
over  the  entire  problem  domain  by  inverting  a  linear  system  formed  in  the  manner 
described  in  the  previous  chapter.  This  linear  system  is  constructed  once  for  each 
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radiative  mechanism  considered  and  for  each  transmission  direction  considered.  For 


instance,  if  three  radiative  mechanisms  and  three  transmissions  directions  were  consider, 
the  linear  system  would  be  constructed  a  total  nine  times.  Note  that  the  coefficients  of 
the  LHS  matrix  and  the  RHS  vector  of  source  terms  do  vary  with  these  different 
realizations  of  the  linear  system  according  to  the  mechanism-specific  emission  and 
absorption  data,  as  well  as  the  different  direction  cosines  formed  between  the  selected 
transmission  directions  and  the  cell  face  within  the  discretized  domain.  The  reader  is 
directed  to  the  provided  pseudo-code  for  details  related  to  the  construction  of  the  linear 
system  described  above. 

While  the  radiative  source  terms  are  calculate  during  the  execution  of  rad  TS  and 
rad_FVM,  they  are  not  in  a  form  which  is  suitable  for  use  in  the  flow  solver.  The  source 
terms  are  thus  made  suitable  for  this  use  by  the  subroutine  rad _couple.  Finally,  these 
source  terms  are  stored  in  the  common  block  rad_common  and  passed  into  the  main 
program  where  they  are  utilized  in  the  subroutine  sourcet.  Since  radiative  intensity  is 
also  a  quantity  of  interest,  it  is  passed  into  the  main  program  and  subsequently  written  to 
an  output  file.  Unlike  the  two  other  sections  of  this  chapter,  which  pertain  to  the  flow 
and  spectroscopic  solvers,  this  section  contains  no  discussion  of  modifications  to  a 
baseline  code.  The  reason  is  that  the  programming  of  the  various  subcomponents  of  the 
radiation  solver  resulted  from  work  conducted  under  the  reported  research  activity.  A 
summary  of  these  programming  activities  is  provided  below  in  Table  4. 
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Table  4.  Summary  of  Code  Development  Activities  pertaining  to  Radiation  Solver. 


Description 

Subroutine 

Appendix  Entries 

Common  Block 

rod common 

Table  14 

Specify  rad  parameters 

rad parameters 

Table  15 

Subroutine  execution 

rod solver 

Table  16 

Tangent  Slab  Solver 

rad TS 

Tables  17-22 

FVM  Solver 

rad FVM 

Tables  23-26 

Banded  LinearSolver 

band 

Table  27 

Calculate  source  terms 

rad_coupied 

Table  28 

Spectroscopic  Solver 

As  challenging  as  the  coupling  of  the  flowfield  and  radiation  solutions  may  be, 
the  tasks  performed  by  the  spectroscopic  solver  SPRADIAN,  as  implemented  in  the 
subroutine  radipac,  are  critical  to  the  accuracy  of  the  results  which  are  obtained  by  the 
overall  solution  method.  The  primary  task  is  the  calculation  of  the  spectral  emission  and 
absorption  coefficients.  This  task  is  supported  by  the  secondary  tasks  of  calculating  the 
state  populations  of  the  atomic  and  diatomic  systems,  together  with  the  calculation  of  the 
transition  probabilities,  line  profiles  and  line  strengths  associated  with  each  transition. 
Figure  5  below  illustrates  the  basic  structure  of  the  radipac. 

The  execution  of  radipac  begins  with  the  passage  of  values  into  the  subroutine^ 
arguments  from  rad_solver.  Setup  of  the  spectroscopic  solver  is  accomplished  with  these 
passed  values  and  the  various  modules  interface jnod,  structure _def_mod,  size_def_mod. 
At  this  point  it,  it  is  worth  noting  that  SPRADIAN  is  a  rather  extensive  code  and 
significant  portions  of  it  are  not  needed  in  the  present  implementation.  Therefore,  only 
those  subroutines  which  have  been  utilized  and  modified  will  be  discussed.  Figure  5 
below  illustrates  the  program  flow  for  the  utilized  components  of  SPRADIAN. 
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Radiation  Solver 
( rad_solver ) 


Figure  5.  Flowchart  for  Spectroscopic  Solver  (SPRADIAN);  bold  indicates  areas  of  the 

code  affected  by  modifications. 

The  subroutine  emis_absb  is  primarily  responsible  for  coordinating  the  execution 
of  the  various  subroutines  which  are  called  in  order  to  calculate  the  emission  and 
absorption  coefficients.  The  first  such  subroutine  is  atomjbb  which  calculates  these 
spectroscopic  coefficients  based  upon  the  calculated  number  density  of  the  internal 
electronic  states  of  atomic  species  and  tabulated  spectroscopic  data,  such  as  transition 
probabilities.  The  second  subroutine  is  diatomjbb ;  it  calculates  the  spectroscopic 
coefficients  based  upon  similar  calculations  but  with  a  few  key  differences.  The  most 
significant  difference  arises  from  the  necessity  of  considering  transitions  between 
electronic,  vibrational  and  rotational  levels  within  the  diatomic  species.  These  transition 
probabilities  are  calculated  based  upon  the  theoretical  developments  presented  in  Chapter 
III  regarding  bound-bound  radiation  in  diatomic  systems. 
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Modifications. 

The  modifications  made  to  the  baseline  spectroscopic  solver,  SPRADIAN,  were 
primarily  concerned  with  passing  data  between  radipac  and  rad_solver  and  with  parsing 
the  total  emission  and  absorption  coefficients  into  their  spectral  components  by 
mechanism.  Listings  of  these  modifications  may  be  found  in  Appendix  A  and  are 
summarized  below  in  Table  5. 


Table  5.  Summary  of  Modifications  Pertaining  to  Spectroscopic  Solver. 


Description 

Subroutine 

Appendix  Entries 

Assign  common  variables 

radipac 

Table  29 

Extract  state  populations 

radipac 

Table  30 

Specify  Tvibs;  added  logic 

emis absb 

Table  31 

Spectral  emis  and  absb 

atom bb 

Table  32 

Spectral  emis  and  absb 

diatom_bb 

Table  33 
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VI.  Results 


The  results  presented  in  this  chapter  are  organized  into  four  sections  consisting  of 
one  section  discussing  the  solution  parameters  as  taken  from  the  conditions  of  the  FIRE  II 
flight  experiment  and  three  sections  corresponding  to  the  phases  of  work  conducted.  The 
first  of  these  three  remaining  sections  details  the  comparison  of  the  multispecies 
multitemperature  and  two-temperature  thermal  model.  The  comparison  was  conducted 
by  examining  the  flow  fields  obtained  by  NH7AIR  and  the  two-temperature  flow  solver 
LAURA.  Both  codes  use  the  same  air  chemistry  model  and  thermophysical  data.  So, 
any  differences  observed  in  both  the  flow  field  and  the  uncoupled  tangent  slab  radiation 
results  are  due  to  the  manner  in  which  internal  energy  is  distributed  among  the  available 
modes. 

The  second  phase  of  work  includes  results  obtained  by  coupling  NH7AIR  and 
SPRADIAN  with  tangent  slab  radiative  transport.  A  comparison  of  these  coupled  results 
with  uncoupled  NH7AIR  results  was  conducted  in  order  to  investigate  the  effects  of 
coupling  the  radiation  source  terms  into  the  flow  solver.  The  effects  on  the  radiation 
solution  observed  are  also  reported. 

The  third  phase  work  corresponds  to  the  development  of  a  finite  volume  method 
for  radiative  transport.  Spectrally  coarse  results  for  three  uncoupled  cases  are  examined 
in  this  final  section, 

and  some  of  the  geometrical  effects  of  the  current  implementation  are  discussed. 
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The  FIRE  II  Flight  Experiment 

Data  from  the  FIRE  II  flight  experiment  have  been  used  to  validate  the  code 
developed  in  the  research  activities  associated  with  this  dissertation.  This  flight 
experiment  was  undertaken  by  NASA  prior  to  the  Apollo  missions  to  investigate  the 
heating  environment  surrounding  vehicles  reentering  the  Earth‘s  atmosphere.  Of 
particular  interest  to  this  test  program  was  the  characterization  of  the  radiance  and  heat 
transfer  rates  on  large-scale  blunt-nosed  bodies.  The  resulting  data  was  intended  for 
comparison  with  ground-based  experiments  and  theoretical  calculations  (Lewis  and 
Scallion,  1965). 

The  spacecraft  configuration  shown  in  Figure  6  included  three  total  radiometers 
(one  on-axis,  one  off-axis  and  one  aft  facing)  as  well  as  a  spectral  radiometer  which  was 
bore  sighted  with  the  on-axis  total  radiometer.  Additionally,  a  calorimeter  monitored  the 
total  heat  load  on  the  forebody.  The  original  data  collection  and  reduction  plan 
anticipated  that  it  would  be  possible  to  determine  the  convective  heat  load  by  subtracting 
the  radiative  heat  flux  measured  by  the  total  radiometer  from  the  total  heat  flux  measured 
by  the  calorimeter.  Figure  7  illustrate  the  various  phases  of  the  flight  experiment. 
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Figure  7.  Mission  Profile  for  the  FIRE  II  Experiment 
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Sequence  of  events  In  Project  Ffrfl  flight  IL 


This  was  all  based  on  the  assumption  that  most  of  the  radiation  emitted  in  the 
shock  layer  would  be  above  the  optical  cutoff  of  the  radiometer  window  X  >  2, 000A . 
However,  it  was  inferred  from  the  subsequent  analysis  of  the  data  collected  by  the  total 
radiometers  and  calorimeter  that  strong  vacuum-ultraviolet  (VUV)  sources  were  present 
in  the  flow  field.  Besides  the  presence  of  strong  VUV  absorption,  it  was  determined  that 
the  primary  source  of  radiative  emissions  were  from  the  near-infrared  lines  of  the  atomic 
flow  species. 

The  reentry  of  the  flight  trajectory  consisted  of  three  distinct  phases,  each  of 
which  corresponds  to  the  ejection  of  one  of  the  layered  heat  shields  depicted  in  Figure  5. 
The  first  phase  occurred  prior  to  the  ejection  of  the  first  heat  shield  beginning  at  a  total 
elapsed  time  of  1631.3  seconds  and  ending  at  1636.5.  During  this  first  period  the  flow 
exhibits  a  range  of  equilibrium  conditions.  The  flow  starts  this  experimental  period  in  a 
state  of  severe  nonequilibrium  and  by  its  conclusion  has  reached  as  state  of  near 
equilibrium.  This  range  of  equilibrium  conditions  makes  this  an  ideal  data  set  for 
validating  a  code  like  the  one  developed  here. 

Solution  Parameters 

The  radiation  along  the  stagnation  line  of  the  1634,  1636,  and  1640.5-second 
trajectory  points  of  the  FIRE  II  experiment  have  been  investigated.  These  trajectory 
points  were  selected  because  of  the  range  of  nonequilibrium  conditions  exhibited:  from 
highly  nonequilibrium  for  the  1634-second  point  to  near  equilibrium  for  the  1640.5- 
second  point.  Table  6  below  contains  the  solution  parameters  at  these  trajectory  points. 
The  freestream  chemical  composition  is  given  in  terms  of  mass  fractions  in  Table  7.  The 
wall  chemistry  was  modeled  using  a  non-catalytic  boundary  condition.  The  grids  used  in 


94 


this  study  both  contained  51x61  nodes  in  the  rotated  plane  of  the  axisymmetric  body. 
The  grid  adaptation  algorithm  of  Gnoffo,  et  al.,  (1993)  was  used  to  place  adequate  points 
in  the  boundary  layer  and  through  the  shock  in  order  to  adequately  resolve  the  gradients 
there.  A  typical  grid  relative  to  the  FIRE  II  vehicle  is  shown  below  in  Figure  8. 


Table  6.  Parameters  for  Flowfield  Solution 


Flow  Field  Parameters 

^elapsed  (^) 

Twall  (K) 

T,nf  (K) 

M 

Re 

Pinf  (Pa) 

1634.0 

615 

195 

40.6 

1.40E+09 

2.08 

1636.0 

810 

210 

38.9 

1.20E+09 

5.16 

1640.5 

1560 

254 

34.4 

8.50E+08 

28.12 

Free  Stream  Mass  Fractions 

^elapsed  (^) 

CN2 

C02 

CN 

co 

CNO,NO+ 

All 

0.767 

0.233 

1.00E-06 

1.00E-06 

1.00E-06 

NFI7AIR  Grid  Dim: 

51x61 

Figure  8.  FIRE  II  Geometry  and  Grid  with  Adaptation  Applied  in  the  Shock  and 

Boundary  Layers 


The  parameters  of  the  spectral  calculations  were  chosen  in  order  to  facilitate  a 
comparison  of  the  different  radiation  solutions  resulting  from  the  multitemperature  and 
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two-temperature  flow  fields.  The  range  of  wavelengths  used  in  this  study  was  2,000- 
40,000  Angstrom  with  100,000  points  used  to  discretize  this  spectral  range.  The 
radiation  bands  and  mechanisms  considered  in  the  present  study  are  listed  in  Table  7. 


Table  7.  Parameters  for  Radiation  Solution 


Species 

Mechanisms 

Key 

State  Transitions 

n2 

Vegard-Kaplan 

VK 

(AX  aaXX) 

n2 

1st  Positive 

1+ 

(b’x^aX) 

n2 

2nd  Positive 

2+ 

(c’n„+A£’nt) 

n2 

Lyman-Birge-Hopfield 

LBH 

(a'ns^xX) 

n2 

Birge-Hopfield  1 

BH1 

(b'n^xX) 

n2 

Birge-Hopfield  2 

BH2 

(b’‘n„^xX) 

02 

Schuman-Runge 

SR 

(bX++xX) 

NO 

P 

Beta 

(a2z+  <r+x2  nr) 

NO 

r 

Gamm 

(b2  nr^x2nr) 

NO 

5 

Delt 

(c2nr^x2n,.) 

NO 

8 

Epsi 

(d2z+  ^x2ur) 

N 

Bound-Bound 

N 

O 

Bound-Bound 

O 

Xmin 

Xmax 

2,000 

40,000 

100,000 
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Comparison  of  Two-Temperature  and  Multitemperature  Models 


Three  trajectory  points  from  the  FIRE  II  experiment  were  chosen  for  use  in  this 
investigation.  This  set  of  trajectory  points  exhibits  a  range  of  nonequilibrium  conditions, 
which  range  from  severe  nonequilibrium  to  near  equilibrium.  The  1634.0-second 
trajectory  point  exhibits  the  highest  degree  of  thermal  nonequilibrium,  while  the  1636- 
and  1640.5-second  trajectory  points  exhibit  progressively  more  equilibrium  behavior.  In 
this  section,  the  effects  of  exchanging  the  two-temperature  model  for  the  multispecies- 
multitemperature  model  are  discussed  according  to  the  characteristic  features  observed  in 
the  flowfield  quantities.  Figure  9  a)  presents  a  comparison  of  the  temperature  profiles 
along  the  stagnation  line  of  the  NH7AIR  and  LAURA  flowfields  at  1634.0  seconds. 
Both  solutions  exhibit  a  shock  stand-off  distance  of  about  7  cm  with  comparable  heavy 
particle  and  electronic  temperatures  in  the  shock  layer. 

However,  some  significant  difference  exists  between  the  two  solutions.  Of  first 
importance,  is  the  fact  that  NH7AIR  predicts  species-vibrational  temperatures  which — 
rising  quickly  within  the  shock — are  far  from  being  at  equilibrium  with  the  electronic 
temperature.  That  vibrational  and  electronic  temperatures  are  in  equilibrium  is  a  key 
assumption  of  the  two-temperature  model.  The  effect  of  allowing  the  species-vibrational 
energy  modes  to  relax  separately  from  the  electronic  modes  is  that  they  are  able  to  do  so 
more  quickly,  according  to  their  relatively  faster  relaxation  times.  This  results  in  a 
predicted  peak  temperature  for  NH7AIR  which  is  about  5,000  K  lower  than  LAURA,  as 
well  as  an  observable  reduction  of  the  shock  thickness.  Another  difference  between  the 
results  of  the  two  solution  methods  becomes  more  noticeable  when  the  remaining  two 
trajectory  points  are  considered.  Figures  9  b)  and  c)  show  that  NH7AIR  consistently 
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predicts  a  reduction  in  peak  temperature  and  shock  thickness.  However,  it  is  apparent 
with  these  two  cases  that  NH7AIR  predicts  a  higher  equilibrium  temperature  in  the  shock 
layer  than  does  LAURA. 


Figure  9.  a)  Comparison  of  Temperature  Profiles:  1634.0  seconds;  s olid  and  dashed  lines 
represent  the  NH7AIR  and  LAURA  data,  respectively. 


Figure  9.  b)  Comparison  of  Temperature  Profiles:  1636.0  seconds 
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Figure  9.  c)  Comparison  of  Temperature  Profiles:  1640.5  seconds 

Figures  10  and  11  compare  the  density  and  pressure  profiles  along  the  stagnation 
streamline  at  each  of  the  three  trajectory  points  considered.  While  both  figures  indicate  a 
substantial  difference  in  shock  strength  and  standoff  distance  at  each  of  the  three  selected 
trajectory  points,  the  solutions  obtained  by  NH7AIR  show  a  reasonable  correlation  in 
pattern  with  those  obtained  by  LAURA  in  terms  of  the  pre-  and  post-shock  flow 
conditions.  The  higher  temperatures  predicted  by  NH7AIR,  together  with  the  post-shock 
pressures  which  are  nearly  identical  to  those  obtained  by  LAURA,  result  in  lower  density 
in  shock-layer  and  thereby  the  greater  shock  standoff  distances  observed  in  the  NH7AIR 
data.  It  is  speculated  that  the  greater  post-shock  temperature  rise  observed  in  the  results 
obtained  under  multispecies,  multitemperature  model  is  the  result  of  a  decreased 
production  of  entropy  relative  to  the  two-temperature  model.  The  transfer  of  energy 
between  energy  modes  which  are  out  of  equilibrium  is  somewhat  analogous  to  the 
transfer  of  energy  to  a  body  from  a  surrounding  heat  bath  (Vincenti  and  Kruger,  1967). 
However,  in  this  case,  energy  in  not  transferred  between  bodies  separated  in  space,  rather 
it  is  transferred  between  energy  modes  separated  by  their  respective  degrees  of  freedom. 
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Here  consider  an  energy  mode  with  a  particular  degree  of  freedom  (e.g.,  simple  harmonic 
oscillator)  which  is  out  of  equilibrium  with  the  surrounding  heat  bath  at  temperature  T  . 
It  can  be  shown  that,  for  an  energy  mode  with  an  energy  content  sufficiently  specified  by 
a  characteristic  temperature  T ,  the  entropy  produced  by  heat  transfer  from  the  heat  bath 
to  this  z-th  energy  mode  is  given  by 


ds<  = 


_L_I 

vTt, 


de. 


(171) 


From  this  expression,  it  is  evident  that  entropy  production  due  to  this  transfer  of  energy 
from  the  heat  bath  to  the  nonequilibrium  energy  mode  is  zero  only  for  the  case  where 
Tt=T  or  det=  0 ,  which  is  the  case  of  thermal  equilibrium.  Therefore,  if  the 

thermodynamic  state  is  closer  to  equilibrium,  these  nonequilibrium  processes  will 
produce  less  entropy;  such  is  the  case  for  the  solutions  obtained  by  the  multispecies- 
multitemperature  thermal  model.  A  reduction  in  entropy  production  means  more  useful 
energy  is  recovered  to  expand  and  heat  the  gas  in  the  shock-layer,  thereby  raising  post¬ 
shock  temperatures,  lowering  density  and  increasing  the  shock  stand-off  distances.  This 
effect  is  more  pronounced  in  the  NH7AIR  results,  which  is  consistent  with  a  reduction  in 
the  production  of  entropy  through  the  shock  wave. 
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Distance  from  Wall  (cm) 

Figure  10.  Comparison  of  Density  Profiles;  solid  and  dashed  lines  represent  NH7 AIR  and 

LAURA  data,  respectively. 


Distance  from  Wall  (cm) 


Figure  11.  Comparison  of  Pressure  Profiles;  solid  and  dashed  lines  represent  NH7 AIR 

and  LA  URA  data,  respectively. 


Figure  12  presents  data  regarding  the  number  densities  calculated  via  NH7AIR 
and  LAURA,  respectively.  Both  depict  physically  realistic  gas  composition  along  the 
streamline,  with  a  high  degree  of  agreement  between  NH7AIR  and  LAURA.  Ahead  of 
the  shock  the  concentrations  reflect  those  specified  at  the  inflow  boundary.  Through  the 
shock,  the  diatoms  N2  and  O2  dissociate  and,  consequently,  the  number  densities  of  these 
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species  drop  several  orders  of  magnitude,  from  about  1020  particles/cm3  down  to  1014 
particles/cm3  at  1634  s  and  1016  particles/cm3  at  1640.5  s.  The  increase  in  the  number 
densities  of  the  atomic  species,  Nn  and  No,  is  the  result  of  strong  dissociation  through  the 
shock,  causing  these  values  to  rise  from  their  freestream  values  to  about  1022 
particles/cm3  in  the  shock-layer.  Higher  atomic  number  densities  occur  at  later  trajectory 
points  due  to  the  higher  freestream  density.  NNo  increases  by  about  five  orders  of 
magnitude  through  the  shock  to  a  typical  value  of  approximately  1017  particles/cm3,  with 
NO  being  present  in  small  numbers  due  to  its  function  as  an  intermediate  reaction 
between  the  diatomic  species  and  the  fully  dissociated  and  ionized  species. 

In  the  post-shock  region  away  from  the  wall,  the  number  densities  of  the  diatomic 
species  either  level  off  or  continue  to  fall.  The  number  density  for  the  only  ionized 
species  considered  in  this  investigation  NNo+  appears  to  rise  quickly  through  the  shock 
and  to  level  off  in  the  shock-layer  around  1019  to  102°  particles/cm3.  NO+  constitutes 
about  0.1%  of  the  flow  in  terms  of  the  total  number  of  particles.  Approaching  the  wall, 
there  is  an  increase  in  the  number  densities  of  the  diatomic  species  then  a  sudden  drop  in 
the  number  densities  of  NO  and  NO+,  corresponding  to  a  sudden  rise  in  the  atomic 
species  at  the  wall.  The  differences  between  the  two  flowfield  solution  methods  in  terms 
of  the  number  densities  of  the  most  prevalent  species  in  the  post-shock  region  (N,  O,  and 
NO+)  are  within  an  order  of  magnitude  or  better  for  each  of  the  selected  trajectory  points. 


102 


Distance  from  Wall  (cm) 

Figure  12.  a)  Comparison  of  N2  Stagnation-Line  Number  Density  Profiles.  Solid  and 
dashed  lines  represent  NH7 AIR  and  LAURA  data,  respectively. 


Distance  from  Wall  (cm) 

Figure  12.  b)  N  Stagnation-Line  Number  Density  Profiles 


0  2  4  6  8 

Distance  from  Wall  (cm) 

Figure  12.  c)  O2  Stagnation-Line  Number  Density  Profiles 


103 


Distance  from  Wall  (cm) 

Figure  12.  d)  O  Stagnation-Line  Number  Density  Profiles 


Distance  from  Wall  (cm) 

Figure  12.  e)  NO  Stagnation-Line  Number  Density  Profiles 


Distance  from  Wall  (cm) 

Figure  12.  f)  NO+  Stagnation- Line  Number  Density  Profiles 


104 


Uncoupled  Radiation. 

In  terms  of  the  radiation  solutions  obtained  from  the  flow  field  results  above,  the 
effects  of  substituting  the  two-temperature  thermal  model  with  the  multispecies- 
multitemperature  model  are  significant.  Stagnation  point  radiative  intensity  at  the 
1634.0-second  trajectory  point  was  estimated  to  be  245.0  (W/cm2-sr)  and  390.0  (W/cm2- 
sr)  for  the  two-temperature  and  multispecies-multitemperature  models,  respectively. 
Variation  of  normal  intensity  is  plotted  along  the  stagnation  line  in  Figure  13;  the  key 
may  be  referenced  to  the  full  mechanism  names  given  previously  in  Table  7.  It  is  readily 
noticeable  that  the  radiation  calculated  from  the  two  flow  field  models  is  very  different 
both  in  terms  of  magnitude  and  spatial  distribution. 
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-  0  -  Epsi 
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Figure  13.  a)  Incoming  Normal  Intensity  along  the  Stagnation  Line:  1634.0  seconds. 


105 


Distance  from  Wall  (cm) 

Figure  13.  b)  Incoming  Normal  Intensity  along  the  Stagnation  Line:  1636  seconds 


Distance  from  Wall  (cm) 

Figure  13.  c)  Incoming  Normal  Intensity  along  the  Stagnation  Line:  1640.5  seconds 


Inspection  of  Figure  13  suggests  that  the  vast  majority  of  the  emission  is  from  the 
bound-bound  transitions  of  the  atomic  species  while  only  a  very  small  fraction  is 
attributable  to  vibrational  bands  of  the  diatomic  species.  It  is  also  noted  that  atomic 
radiative  emission  occurs  at  a  higher  rate  in  the  shock  for  the  NH7AIR  data.  This  effect 
is  due  to  the  higher  dissociation  rates  there,  resulting  from  the  much  higher  vibrational 
temperatures.  Table  8  presents  a  comparison  of  the  radiative  intensity  reaching  the  wall 
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along  the  stagnation  line  for  the  uncoupled  cases  considered  above.  The  total  radiative 
intensities  presented  here  have  been  integrated  over  a  spectral  range  of  0. 2-4.0  pm  in 
order  to  compare  them  with  the  total  radiometer  measurements  collected  aboard  the  FIRE 
II  flight  experiment.  The  reader  is  reminded  that  the  cases  considered  here  are  uncoupled 
from  the  radiative  heat  transfer  mechanisms  and  thus  do  not  take  into  account  the  very 
significant  effect  of  radiation  cooling  which  can  drastically  lower  the  radiative  intensity 
within  the  flowfield.  As  will  be  shown  in  the  next  section,  the  effect  of  radiation  cooling 
can  reduce  these  values  by  an  order  of  magnitude  or  more  at  the  flow  conditions. 


Table  8.  Radiative  Intensity  at  the  Stagnation  Point  for  Uncoupled  Cases  with  Flight  Data 
from  the  FIRE  II  Experiment.  Results  integrated  over  a  spectral  range  of  0. 2-4.0  pm. 


Time 

(s) 

NH7AIR 

(W/cm2-sr) 

LAURA 

(W/cm2-sr) 

FIRE  II 
(W/cm2-sr) 

1634.0 

68.0 

50.1 

1.3 

1636.0 

180.0 

114.0 

5.0 

1640.5 

390.0 

245.0 

35.0 

Comparison  of  Uncoupled  and  Coupled  Radiative  Transport  Results 

The  observations  made  above  regarding  the  comparison  of  uncoupled  radiative 
transport  solutions  obtained  under  the  two-temperature  and  multispecies- 
multitemperature  thermal  models  reveal  that  although  solutions  under  the  multispecies- 
multitemperature  model  may  bear  a  certain  resemblance  to  their  two -temperature 
counterparts  it  terms  of  peak  temperature,  shock  stand-off  distance  and  chemical 
composition,  marked  differences  were  easily  distinguished  in  terms  of  the  nonequilibrium 
distribution  of  energy  among  the  various  energy  modes.  The  most  noticeable  difference 
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was  the  apparent  violation  of  the  assumption  inherent  in  the  two -temperature  model, 
namely,  that  the  vibrational  and  electronic  energy  manifolds  are  far  from  equilibrium 
within  both  the  shock  and  relaxation  zone.  This  effect  is  perhaps  not  surprising, 
considering  that,  due  to  the  disparity  in  mass  between  electrons  and  heavy  particles,  the 
energy  exchanges  which  occur  between  the  two  are  relatively  inefficient  as  compared  to 
energy  exchanges  between  heavy  particles  and  diatomic  molecules  (Park,  1991).  In 
Figure  14  a)  -  c),  the  disparity  between  vibrational  and  electronic  temperatures  is 
observed  quite  readily.  The  vibrational  temperatures  rise  throughout  the  very  diffuse 
shock  and  equilibrate  with  the  heavy  particle  temperature  downstream  of  the  shock,  while 
the  electron  temperature  climbs  slowly  through  both  the  shock  and  subsequent 
downstream  region,  finally  equilibrating  with  the  heavy  particle  and  vibrational 
temperatures  just  before  reaching  the  wall. 


Figure  14.  Stagnation-Line  Temperatures  from  Coupled  NH7AIR-SPRADIAN  Solutions: 

a)  1634.0  seconds 
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Leaving  aside  the  results  of  the  uncoupled  investigation,  the  coupled  solutions  are 
now  examined  as  obtained  by  the  multispecies-multitemperature  nonequilibrium  flow 
solver.  First  of  all,  certain  features  of  the  solution  were  striking.  One  might  expect  a 
drop  in  all  temperatures  in  the  shock  layer,  a  higher  rate  of  recombination  in  the  shock 
layer,  and  a  reduction  in  shock  standoff  distance.  Instead,  results  show  a  dramatic  change 
in  the  nonequilibrium  energy  distribution  within  the  flowfield  and  almost  no  change  in 
flow  composition.  The  most  significant  effect  of  coupling  the  radiation  source  terms  with 
the  nonequilibrium  flow  solver  was  the  dramatic  drop  in  electron  temperature.  This  drop 
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in  electron  temperature  was  accompanied  by  a  modest  rise  in  heavy  and  species- 
vibrational  temperatures. 

While  coupling  radiation  source  terms  into  the  flow  solver  resulted  in  a  significant 
reduction  of  the  electronic  temperature  within  the  flow  field,  this  had  a  negligible  effect 
on  flowfield  composition  throughout  most  of  the  solution  domain.  Notwithstanding  this 
result,  some  small  variations  in  flow  composition  were  noted  for  a  few  of  the  coupled 
cases  in  the  regions  of  the  flow  near  the  wall  and  traversing  the  shock.  For  instance, 
Figure  15  shows  that  the  coupled  result,  for  the  1634.0-second  trajectory  point,  exhibits  a 
faster  ionization  rate  traversing  the  shock  due  to  the  higher  electronic  temperature  there. 


Figure  15.  Stagnation-Line  Number  Density  Profiles  of  NO+:  1634.0-second;  solid  and 
dashed  lines  represent  the  coupled  and  uncoupled  NH7AIR-SPRADIAN  cases, 

respectively. 

Also,  it  is  noted  that  coupling  seems  to  have  had  the  effect  of  slowing  the 
recombination  of  N  and  O  to  produce  O2  and  NO  near  the  wall,  for  the  1636-second 
coupled  case,  as  shown  in  Figures  16  a)  and  b),  consistent  with  the  higher  electronic 
temperature  observed  in  the  uncoupled  cases. 
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Figure  16.  a)  Stagnation-Line  Number  Density  Profiles  of  CL:  1636  seconds.  Solid  and 
dashed  lines  represent  the  coupled  and  uncoupled  NH7AIR-SPRADIAN  cases, 

respectively. 


Figure  16.  b)  Stagnation-Line  Number  Density  Profiles  of  NO:  1636  seconds 
Coupled  Radiation. 

Results  presented  here  are  delineated  according  to  two  parameters.  The  first 
parameter  is  the  total  integrated  radiative  intensity,  calculated  using  the  tangent-slab 
approximation  to  radiative  transport  equation.  Figures  17  a)  -  c)  shows  the  profiles  of 
integrated  intensity  for  the  selected  trajectory  points.  Intensity  profiles  have  been  plotted 
on  a  logarithmic  scale  in  order  to  illustrate  the  diverse  range  of  contributions  from  the 
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various  participating  radiative  mechanisms  to  the  total  radiative  intensity.  This  figure, 
which  is  analogous  to  the  uncoupled  results  in  Figure  13,  readily  shows  that  the 
dominating  radiative  mechanisms  are  from  the  line  emissions  of  the  atomic  species,  most 
notably  nitrogen,  which  generally  accounts  for  as  much  as  90%  of  the  total  radiation  in 
the  cases  investigated.  The  molecular  band  mechanisms  contribute  much  less  to  the  total 
radiation  relative  to  the  atomic  line  radiation,  due  largely  to  being  much  fewer  in  number 
relative  to  atomic  species.  This  point  is  illustrated  by  the  rise  in  radiative  intensity 
approaching  the  wall  where  recombination,  together  with  a  sufficiently  high  electronic 
temperature,  affects  a  marked  rise  in  net  radiative  emission  from  the  molecular 
mechanisms. 


i°T 


0  1  2  3  4  5  6  7 

Distance  from  Wall  (cm) 


Figure  17.  a)  Incoming  Normal  Intensity  along  the  Stagnation  Line:  1634.0  seconds;  solid 
and  dashed  lines  represent  the  coupled  and  uncoupled  NH7AIR-SPRADIAN  cases, 

respectively. 
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In  general,  radiative  emission  is  a  strong  function  of  upper  electronic  state 
populations;  as  such,  these  populations  serve  as  the  second  parameter  along  which  the 
radiative  results  may  be  examined.  Figures  18  a)  and  b)  show  a  representative  pair  (i.e., 
coupled  and  uncoupled  solution)  of  state  populations  for  N. 
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Distance  from  Wall  (cm) 

Figure  18.  a)  Nonequilibrium  Group  Populations  ofN:  Uncoupled  NH7 AIR,  1634.0 

seconds 


012345678 
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Figure  18.  b)  Coupled  NH7AIR-SPRADIAN,  1634  seconds 


Two  trends  are  evident  looking  at  the  differences  in  the  state  populations  of  the 
various  species  at  each  of  the  cases.  First,  the  populations  of  the  ground  states  are  fairly 
similar  (i.e.,  same  order  of  magnitude)  between  coupled  and  uncoupled  cases.  Second, 
the  state  populations  of  the  upper  energy  states  show  orders  of  magnitude  differences 
between  coupled  and  uncoupled  solutions.  This  effect  is  due  to  the  change  in  the  electron 
temperature.  Coupled  solutions  exhibit  a  radiative  cooling  effect  of  the  electrons  which 
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redistributes  the  internal  energy  to  the  lower  energy  levels  for  radiating  species.  Plots  of 
the  state  populations  for  the  remaining  species  and  trajectory  points  are  listed  in 
Appendix  B. 

Table  9  summarizes  these  results  by  way  of  comparison  with  experimental  data 
obtained  from  the  FIRE  II  flight  experiment.  As  expected,  coupling  the  radiation  and 
flowfield  solutions  had  the  effect  of  significantly  reducing  the  amount  of  radiative  energy 
present  in  the  solution  domain,  thus  reducing  the  amount  of  radiation  incident  on  the 
wall.  Given  the  uncertainties  involved  with  both  the  flight  data  collected  from  the  FIRE 
II  flight  experiment  and  the  thermophysical  data  available  for  these  sorts  of 
computations,  the  coupled  NH7AIR-SPRADIAN  results  agree  quite  well  with  the 
experimental  data.  The  best  agreement  was  observed  in  the  near  equilibrium  conditions 
of  the  1640.5-second  trajectory  point,  and  the  least  agreement  was  observed  in  the  severe 
nonequilibrium  conditions  of  the  1634.0-second  trajectory  point.  This  trend  in  the  errors 
is  typical  of  the  results  obtained  by  other  researchers  (Johnston,  2006).  However,  it  is 
key  to  note  that  the  results  presented  here  not  only  agree  reasonably  well  the  FIRE  II  data 
but  underpredict  the  amount  of  radiation  observed.  The  significance  of  this  result  arises 
in  light  of  the  fact  that  the  radiometer  windows  on  the  FIRE  II  vehicle  were  recessed  into 
the  heat  shields,  thereby  capturing  some  radiating  ablation  products.  This  trapped 
ablation  material  contaminated  the  radiometer  data  to  some  unknown  degree,  resulting  in 
reporting  of  inflated  intensity  measurements  (Greendyke,  2011).  Logically,  the 
contribution  of  the  air  species  in  the  flow  field  to  the  radiative  intensity  must  be  lower 
than  the  reported  values.  Thus,  the  agreement  between  the  coupled  NH7AIR-SPRADIAN 
results  reported  here  and  the  FIRE  II  data  is  better  than  it  may  initially  seem,  especially 
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for  the  case  of  severe  nonequilibrium  where  calculated  intensities  tend  to  be  higher  than 
those  reported  in  the  FIRE  II  data  or  from  comparable  results  in  the  literature  as  reported 
by  Johnston  (2006). 

Table  9.  Radiative  Intensity  at  the  Stagnation  Point  for  Coupled  NH7AIR-SPRADIAN 
FIRE  II  Cases.  Results  integrated  over  a  spectral  range  of  0. 2-4.0  pm. 


Time 

(s) 

Uncoupled 

(W/cm2-sr) 

Coupled 

(W/cm2-sr) 

Exp. 

(W/cm2-sr) 

Literature  (lo-hi) 
(W/cm2-sr) 

1634.0 

68.0 

0.5 

1.3 

03-2.6 

1636.0 

180.0 

2.8 

5.0 

2. 6-6. 8 

1640.5 

390.0 

30.0 

35.0 

12.0-39.0 

Finite  Volume  Method  for  Radiative  Transport 

The  finite  volume  method  for  radiative  transport  (FVMR)  within  the  flow  field 
yields  a  full  3 -dimensional  solution  to  the  radiative  transport  equation  which  conserves 
radiative  energy.  This  property  of  conserving  radiative  energy  makes  the  FVMR  a 
desirable  method  for  calculating  source  terms  for  coupling  with  a  flow  solver — an 
especially  important  consideration  for  strongly  absorbing  media  such  as  air  in  the  VUV 
spectra.  However,  without  a  parallelized  solution  algorithm,  obtaining  solutions  using 
this  method  is  computationally  prohibitive.  Since  it  was  outside  of  the  scope  of  this  work 
to  parallelize  the  computer  code,  only  approximate  results  using  the  FVMR  approach  are 
presented  here.  These  results  are  approximated  by  considering  the  radiative  transport 
resulting  from  emission  and  absorption  coefficients  which  have  been  calculated  with  a 
very  coarse  spectral  resolution.  A  spectral  resolution  of  10,000  grid  points  has  been  used 
here  in  contrast  to  the  100,000  grid  points  used  in  the  preceding  sections.  Also,  the 
solutions  presented  here  were  calculated  using  the  thermodynamic  variables  obtained 
from  the  flow  fields  in  the  previous  section.  Trying  to  calculate  coupled  solutions  with 
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such  a  coarse  spectral  resolution  would  not  have  yielded  any  result  more  meaningful  than 
those  presented  here.  Finally,  the  radiative  transport  equation  was  solved  using  the 
spectrally  coarse  coefficients. 

Table  10  summarizes  the  results  obtained  via  the  coarse  FVMR  calculation 
described  above.  These  FVMR  results  are  presented  alongside  tangent  slab  results  at  the 
same  level  of  spectral  resolution  for  comparison.  There  appears  to  be  fairly  good 
agreement  between  the  two  methods  for  the  1634.0  and  1636.0  second  trajectory  points. 
It  is  suspected  that  the  unusually  high  value  for  the  FVMR  at  1640.5  seconds  is  due  to 
geometrical  effects. 

Table  10.  Integrated  Radiative  Intensity  at  the  Stagnation  Point  for  Uncoupled  Cases 
(NH7AIR-SPRADIAN).  Results  integrated  over  a  spectral  range  of  0. 2-4.0  pm  with 

N=10,000  spectral  grid  points. 


Time 

Tangent  Slab 

FVMR 

(s) 

(W/cm2-sr) 

(W/cm2-sr) 

1634.0 

2.2 

2.3 

1636.0 

21.0 

24.6 

1640.5 

187.6 

671.0 

The  reader  will  recall  that  6  transmission  directions  were  considered  in  this  work. 
The  first  two  plots  presented  in  Figures  19  a)  and  b)  are  for  the  positive  and  negative  z 
direction  transmission  directions,  respectively.  The  z-axis  runs  along  the  line  of 
symmetry  and  is  positive  in  the  direction  away  from  the  body.  These  solutions  are 
roughly  analogous  to  the  tangent  slab  solution  in  the  stagnation  region.  The  second  set  of 
plots  presented  in  Figures  19  c)  and  d)  are  for  the  positive  and  negative  r  direction 
transmission  directions.  The  /"-axis  runs  radially  from  the  line  of  symmetry  out  to  the 
farfield.  These  two  solutions  exhibit  the  geometric  effects  discussed  in  the  previous 
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chapter  related  to  the  calculation  of  radiative  transmission  via  the  FVMR  within  an 
axisymmetric  wedge.  Finally  the  solution  in  the  6  direction  is  presented  in  Figure  18  e). 
The  6*-axis  is  out  of  the  plane  in  the  figures  below  and  has  a  circumferential  orientation  in 
the  coordinate  system.  The  reader  will  recall  that  it  is  not  possible  to  use  the  FVMR  to 
calculate  a  solution  for  transmission  directions  which  include  a  component  in  the  0 
direction  with  the  grid  topology  used  here.  In  order  to  calculate  the  radiation  transmitted 
in  these  directions  it  is  necessary  to  use  the  complete,  3 -dimensional  domain.  Therefore, 
the  result  presented  below  simply  represents  the  net  intensity  radiated  from  the  wedge 
under  consideration. 
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Figure  19.  b)  Spectrally  Integrated  z- Direction  Intensity,  V  ,  from  Uncoupled  FVMR 

Solution,  1634  seconds 


Figure  19.  c)  Spectrally  Integrated  r+  Direction  Intensity,  /'  + ,  from  Uncoupled  FVMR 

Solution,  1634  seconds 
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Figure  19.  e)  Spectrally  Integrated  6  Direction  Intensity,  1 9 ,  from  Uncoupled  FVMR 

Solution,  1634  seconds 
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V.  Conclusions 


As  the  United  States,  its  allies  and  its  foes  continue  to  pursue  the  development  of 
new  hypersonic  systems,  the  computational  modeling  of  phenomena  associated  with 
hypersonic  flight  will  play  a  key  role  in  unlocking  the  physical  understanding  requisite  to 
their  design,  manufacture  and  deployment.  Furthermore,  given  the  development  and 
weaponization  of  such  systems,  highly  accurate  modeling  of  radiating  shock  layers  may 
provide  the  critical  MASINT  data  which  will  enable  the  timely  detection  and 
neutralization  of  threats  of  this  kind  to  the  US  and  its  allies. 

Radiation  modeling  has  been  extensively  studied,  particularly  with  respect  to  the 
atmospheric  reentry  of  spacecraft,  as  exemplified  by  the  breadth  and  depth  of  literature 
on  the  subject.  Numerous  computer  codes  have  been  developed  for  modeling  the 
radiation  produced  in  these  situations.  The  level  of  approximation  accepted  in  these 
computer  codes  has  varied  from  those  utilizing  simple  band  models  in  order  to 
characterize  the  spectral  variation  of  the  radiative  transport  properties  to  highly 
sophisticated,  computationally  expensive  line-by-line  methods.  In  the  past  couple  of 
decades  especially,  what  all  these  methods  have  shared  in  common  has  been  the 
utilization  of  the  two-temperature  model  of  thermal  nonequilibrium.  The  present  work 
has  sought  to  advance  the  state-of-the-art  by  proposing  a  more  detailed  model  of 
nonequilibrium,  namely  the  multispecies,  multitemperature  model.  In  this  dissertation,  a 
complete  computational  method  has  been  developed  around  the  line-by-line  radiation 
solver  SPRADIAN  and  the  sophisticated  nonequilibrium  flow  solver  NH7AIR  which 
implements  this  multispecies,  multitemperature  model. 
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The  first  phase  of  code  development  utilized  the  standard  tangent  slab  method  of 
solving  the  radiative  transport  equation.  Results  were  obtained  first  without  coupling 
radiative  effects.  These  uncoupled  results,  obtained  utilizing  the  multitemperature  flow 
solver  NH7AIR,  were  compared  to  the  uncoupled  results  obtained  utilizing  the  two- 
temperature  flow  solver  LAURA.  The  result  of  this  comparison  was  to  show  that,  in  the 
flow  fields  of  the  FIRE  II  cases  which  were  examined,  the  two-temperature  model  does 
not  describe  the  nonequilibrium  processes  involved  as  well  as  the  multitemperature 
model.  The  two-temperature  model  accounts  for  the  redistribution  of  internal  energy 
among  the  vibrational,  electronic  and  free  electron  manifolds,  but  not  with  as  much 
fidelity  as  the  multi-species,  multi-temperature  model.  This  lack  of  fidelity  in  previous 
methods  has  significant  implications  for  the  characterization  of  the  spectral  features  of 
radiating  gases  modeled  in  reentry  shock  layers,  since  the  radiative  properties  of  the  flow 
field  depend  in  a  strongly  nonlinear  fashion  upon  the  temperatures  which  are  calculated 
as  a  result  of  these  nonequilibrium  models. 

Next,  the  tangent  slab  method  was  implemented  within  the  flow  field-radiation 
solver  in  a  coupled  manner  and  validated  against  data  collected  during  the  FIRE  II  flight 
experiment.  The  coupled  implementation  of  the  NH7AIR  and  SPRADIAN  with  the 
tangent  slab  method  dramatically  illustrated  the  effects  of  radiative  cooling  in  the 
modeling  of  reentry  shock  layers.  Furthermore,  excellent  agreement  was  obtained  with 
the  FIRE  II  experimental  data,  especially  for  the  severe  nonequilibrium  conditions  of  the 
1634.0-second  trajectory  point.  An  1 1 -species  air  chemistry  model  would  likely  improve 
these  results  still  further. 
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The  second  phase  of  code  development  focused  on  the  implementation  of  a 
suitable  FVM  scheme  for  solving  the  radiative  transport  equation.  This  FVMR  solution 
method  was  successfully  developed  and  implemented  in  an  uncoupled  fashion  within  the 
developed  computer  code.  The  results  compared  with  those  obtained  from  the  tangent 
slab  method.  The  FVMR  of  calculating  radiative  intensity  is  extremely  memory 
intensive  because  of  the  extensive  linear  system  created  by  attempting  to  resolve  both  the 
spatial,  directional  and  spectral  contributions  of  the  radiation  solution.  It  is  necessary  to 
parallelize  the  FVMR  in  order  to  use  the  level  of  spectral  resolution  needed  in  order  to 
calculate  an  accurate  coupled  flow  field-radiation  solution.  In  an  effort  to  present  some 
manner  of  result,  spectrally  coarse,  uncoupled  calculations  were  performed  on  the  FIRE 
II  flow  fields  in  order  to  obtain  both  tangent  slab  and  FVMR  solutions.  Comparing  these 
two  solutions  yielded  a  reasonable  amount  of  agreement.  However,  in  the  course  of 
analyzing  the  results  of  the  FVMR  scheme  some  undesirable  geometric  effects  were 
observed  which  indicate  that  this  method  would  be  more  appropriately  applied  in  a  fully 
three-dimensional  radiation  grid  rather  than  the  axisymmetric  wedge  used  for  the  flow 
solver. 

The  further  development  of  the  FVMR  should  continue.  Since  it  is  based  on  a 
conservation  law,  it  should  yield  more  physical  results  than  the  tangent  slab  method  when 
coupled  with  a  flow  solver.  This  improvement,  together  with  the  improvements  afforded 
by  the  multispecies  multitemperature  thermal  model,  will  ultimately  result  in  superior 
coupled  flowfield-radiation  solutions  compared  with  present  capabilities.  This  improved 
modeling  capability  may  one  day  aid  in  the  development  of  a  high-performance  reusable 
space  access  vehicles  or  hypersonic  cruise  missile  technology.  Alternately,  they  may 
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serve  to  populate  a  database  of  MASINT  signatures  used  to  identify  incoming  threats.  In 
this  case,  the  accuracy  of  these  methods  could  mean  the  difference  between  a  castrophic 
surprise  attack  by  a  stealthy  hypersonic  weapon  system  and  a  successful  defense  against 
such  threats. 

Recommendations. 

In  order  to  further  improve  the  solution  method  presented  in  this  dissertation,  a 
few  key  recommendations  are  here  made  for  the  consideration  of  those  who  may  desire  to 
develop  this  method  further.  The  first  recommendation  is  that  the  work  of  updating  the 
chemistry  model  be  undertaken.  Updating  the  chemistry  model  to  an  11 -species  air 
model  will  enable  a  more  accurate  calculation  of  the  flow  composition  at  the  conditions 
of  interest.  The  additional  ionization  processes  will  also  have  the  effect  of  lowering  the 
post-shock  temperatures,  which  in  turn  will  have  an  effect  on  the  amount  of  radiation 
produced  and  coupled  into  the  flow  field  solution. 

The  second  recommendation  is  that  a  careful  study  of  the  combined  flow  field  - 
radiation  solver  be  undertaken  in  order  to  determine  optimum  method  by  which  to 
parallelize  the  code.  There  are  many  time-intensive  calculations  within  the  radiation 
solver  which  are  physically  independent  and  would  lend  themselves  well  to  a  parallelized 
implementation. 

Finally,  it  is  recommended  that  work  be  done  in  order  to  develop  a  method  by 
which  the  thermodynamic  flow  quantities  may  be  interpolated  onto  a  separate  grid  which 
has  been  optimized  for  solving  the  radiative  transport  equations.  After  obtaining  the 
radiation  solution  on  this  optimized  grid,  the  source  terms  could  then  be  interpolated  back 
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onto  the  flow  field  grid  and  in  this  way  be  coupled  into  the  flow  solver.  Additionally, 
attention  should  also  be  given  to  enhancing  the  stability  of  the  radiative  coupling. 
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Appendix  A:  Selected  Listings  of  Computer  Code 


Table  11.  Parameters  Directing  Radiation  Solution  added  to  Read  Statements 


Subroutine  dating  bold  text  indicates  modifications  to  baseline  code. 


OPEN (UNIT=XX, FILE=INPUT) 

READ (XX, * )  i-index,  j -index 
READ (XX, * )  flow  solver  parameters 
READ (XX, * )  grid  adaptation  parameters 
READ (XX, * )  flow  field  reference  values 

READ (XX,*)  RADINT ,RADRE AD , IMETHOD , I S TAG , RADITV, RADOUT 

READ (XX,*)  i/o  parameters 
CLOSE (XX) 


RADINT  =  X 
RADREAD  =  X 
IMETHOD  =  X 
I STAG  =  X 
RADITV  =  X 
RADOUT  =  X 


!  1  =  radiation  solver  on 

!  1  =  read  prev  rad  soln 

!0  =  Tangent  Slab  Solver;  0  =  FVMR 

!  1  =  Perform  stagnation  line  calculations 

!N  =  Num  iters  b/w  calls  to  rad_solver 

! 1  =  output  rad  solution  for  restart 


Table  12.  Call  to  Radiation  Solver  within  Main  Loop 


Program  main;  bold  text  indicates  modifications  to  baseline  code;  pseudo-code. 

!  MAIN  LOOP 

DO  N=NSTART , NEND 

CALL  subroutine_l 
CALL  subroutine_2 
CALL  subroutine_3 
NTIME=NTIME+DT 
CALL  subroutine_4 
IF (RADINT . eq . 1 ) THEN 

IF (MOD (N,INT (RADITV) ) .eq.O  .and.  (N . NE . NEND) ) THEN 
write ( * , * )  ' call  rad_solver ' 

CALL  RAD_S OLVER (limits , X , Y , T , T_vib_s , rho_i , c_i , parameters ) 

END  IF 

END  IF 

IF ( ( INT_5 . eq . 1 ) .and. ( INT_7a . eq . 1 ) ) THEN 
IF (INT_8 .eq. 0) THEN 
CALL  subroutine_5 
ENDIF 

ENDIF 
END  DO 
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Table  13.  Calculation  of  Source  Terms 


Subroutine  source;  pseudo-code 


*****source  terms  due  to  radiation 
if (radint . eq . 1 ) then 
du_ev_02  =  du_ev_02-DTDG*Q_rad_02 
du_ev~N2  =  du_ev_N2-DTDG*Q_rad_N2 
du_ev_NO  =  du_ev_NO-DTDG*Q_rad_NO 
du_eel  =  du_eel  -DTDG*Q_rad_el 
du_tot  =  du_tot  -DTDG*Q_rad_tot 
end  if 


Table  14.  Common  Block  Used  by  main  and  radjsolver 


Common  Block  mdjwmon 


j *****variable  belonging  to  common  block  /rad/******************************* 
integer  : :  imech,  nmech 
character (4)  ::  mech_name ( 68 ) 

real*8,  pointer  : :  spect_emis spect_absb alpha_vib ( : , : ) 
real*8,  pointer  ::  wave_length ( : ) 

common/ rad/ spect_emis ,  spect_absb, wave_length, nmech,mech_name, alpha_vib 

i **************************************************************************** 


! *****variable  belonging  to  common  block  /radsoln/******************* 


real*8 , 

pointer  : 

tot  emis 

tot  absb 

real*8 , 

pointer  : 

spect  emis  m ( : , 

spect  absb  m ( : 

real*8, 

pointer  : 

wave  length  m(: 

: ,  : ) , alpha  vib  m ( : ,  : 

real*8 , 

pointer  : 

:  Q  rad  tot 

) 

real*8 , 

pointer  : 

norm  int 

:) 

real*8 , 

pointer  : 

Q_rad_s 

real*8, 

pointer  : 

Q  rad  vib 

l  Ni_02  (:,:,:)  ,  Ni_N2  ( 

H- 

1 

o 

1 

-H 

, ,  Ni_NO ( : 

common/ radsoln/ norm_int , Q_rad_tot , Q_rad_s , Q_rad_vib,  tot_emis , tot_absb, 
1  Ni_02 ,  Ni_N2,  Ni_NO,  Ni_N,  Ni_0 

i  ********************************************************************** 


)  , 


Table  15.  Calculate  nmech  and  Store  Values  for  mech  name 


Module  rad parameters _ 

nmech  =  1 

mech_name (1)  =  'Tot.' 
do  i=l , 18 

if (atom_rads (1, i) . ne . '  ' . and . atom_rads (2 , i ) . eq . ' bb ' )  then 

nmech  =  nmech  +  1 

mech_name (nmech)  =  atom_rads ( 1 , i ) 
end  if 

end  do 
do  i=l ,  40 

if (diatom_bands ( 1 , i ) . ne . '  ')  then 

nmech  =  nmech  +  1 

mech_name (nmech)  =  diatom_bands (2 , i ) 
end  if 

end  do 
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Table  16.  Control  Sequence  for  Executing  Subroutines  Local  to  rad_solver 


Subroutine  rad_solver 


!  CALL  RTE  SOLUTION  ALGORITHM 

if (METHODINT.eq. 1) then 
!  write (*,*)  'called  rad_TS ' 

call  rad_TS (STAGINT, ILD, JLD,X, Y, T, TvN2 , Tv02 , TvNO, TvNOpl, Tel,  & 

&  numN2,  num02 ,  numN,  numO, numNO, numNOpl ,  numel , numatom, nummol , numhvy, molwt ) 
else 

!  write  (*,*)  'called  rad_FVM' 

!  call  rad_FVM (STAGINT, ILD, JLD,X,Y,T, TvN2 , Tv02 , TvNO, TvNOpl, Tel, 

& 

!  &  numN2,  num02 ,  numN,  numO, numNO, numNOpl ,  numel , numatom, nummol , numhvy, molwt) 

end  if 

I  •k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

!  CALL  SOURCE  TERM  ALGORITHM 

call  rad_couple (ILD, JLD)  (See  Table  27) 

write (*,*)  'exiting  rad_solver' 

RETURN 

CONTAINS 


Table  17.  Basic  Outline  of  Tangent  Slab  Radiation  Solver;  important  aspects  of  the 
subroutine  radTS  are  further  described  in  the  tables  indicated  below. 


Subroutine  radTS _ 

subroutine  rad_TS 

-Declare  Variables 
-Initialize  Variables 

do  i=l , ILD 
do  dir=l,2 

spect_int  =  0.0 
spect_int_old  =  0.0 


do  j =a, b, increment 

-Calculate  emis  and  absb  (See  Table  18) 
do  imech=l , nmech 

-Calculate  intensity  (See  Table  19) 

-Calculate  source  terms  for  idir  (See  Table  22) 
end  do  lover  imech 

end  if 

-Calculate  total  intensity  (See  Table  21) 
end  do  lover  j 

-Calculate  total  source  terms  (See  Table  22) 
end  do  lover  dir 
END  DO 

end  subroutine  rad  TS 
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Table  18.  Calculate  Emission  and  Absorption  Coefficients 


Subroutine  rad_TS 

if (dir . eg . 1 ) then 

call  radipac ( 0 

. 0 , nnode, 0.0, method,  0.0, 0.0, 0.0,  & 

& 

wavmin , wavmax , 

nwave, avg  num, & 

& 

atom  noneqs,  atom  rads,  diatom  noneqs,  diatom  bands , triatom  bands, & 

& 

T(i,  j)  ,T(i,  j)  , 

TvN2  (i, 

j ) , Tel ( i , j ) , Tv02 ( i , j ) , TvN2 ( i , j ) , TvNO ( i , j ) , & 

& 

0.0, 0.0, 0.0,0. 

0,0. 0,0 

.0,0.0,0.0,0.0,0.0,0.0,& 

& 

numN ( i , j ) , numN2 ( i , j ) , 

numN2p ( i , j ) , numel ( i , j ) , numNO ( i , j ) , numNp ( i , j ) , & 

& 

numO ( i , j ) , num02 ( i , j ) , 

0.0, numOp ( i , j ) , numhvy ( i , j ) , numatom ( i , j ) , nummol ( i , j ) ,  & 

& 

molwt (i, j ) , Ni 

02  < : , i. 

j ) , Ni_N2 ( : , i , j ) , Ni_NO ( : , i , j ) , Ni_N ( : , i , j ) , Ni_0 ( : , i , j ) ) ! , 

spect  emis  m ( : 

, :,i, j) 

=  spect  emis 

spect  absb  m ( : 

,  :/i,  j) 

=  spect  absb 

wave  length  m( 

j) 

=  wave  length 

alpha  vib  m ( : , 

=  alpha  vib 

end  if 

spect  emis  m ( 1 , : 

,i,j)  = 

0.0D0 

spect  absb  m ( 1 , : 

,i,j)  = 

0.0D0 

do  imech=2 , nmech 

spect  emis  m ( 1 , : 

/ if  j  )  = 

:  spect  emis  m(l,:,i,j)  +  spect  emis  m (imech, :, i, j ) 

spect  absb  m ( 1 , : 

,i,j)  = 

:  spect  absb  m(l,:,i,j)  +  spect  absb  m (imech, :, i, j ) 

end  do 

if (Tel (i, j )  .It. 

2000 . 0D0 ) then 

! spect  emis  m ( : , 

hi,  j) 

=  0.0D0 

spect  absb  m ( : , : 

,i,j)  = 

0.0D0 

end  if 

Table  19.  Calculate  Spectral  and  Normal  Intensities 


Subroutine  rad_TS _ 

spect_int_old=spect_int 

do  imech  =  1,  nmech 
do  m  =  1 , nwave 

if (depth*spect_absb_m (imech,m, i, j ) . gt . 1 . Oe-4 )  then 
blam=spect_emis_m (imech, m, i, j ) / spect_absb_m (imech, m, i, j ) 
spect_int (imech, m, dir) =blam- (blam-spect_int_old (l,m, dir) ) *& 

&  exp (-spect_absb_m (imech, m, i, j ) *depth) 

else 

spect_absb_m (imech, m, i, j )  =  0.0 
spect_int (imech, m, dir )  =  spect_int_old ( 1 , m, dir )  + 
spect_emis_m (imech, m, i, j ) * depth 
end  if 

spect_int (imech, m, dir) =spect_int (imech, m, dir) -spect_int_old (l,m, dir) 
spect_int (imech, m, dir) =spect_int_old (imech, m, dir) +spect_int (imech, m, dir) 
if (spect_emis_m (imech, m, i, j ) .lt.l.0E-20)  spect_emis_m (imech, m, i, j )  = 

0.0D0 

end  do 

! Calculate  Normal  Intensities 
call 

simpson (norm_int (imech, i, j , dir) , wave_length_m ( : , i, j ) , spect_int (imech, : , dir) , nwave, ier) 
norm_int (imech, i, j , dir)  =  norm_int (imech, i, j , dir) *1 . Oe-4 
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Table  20.  Calculate  Partial  Source  Terms;  Calculate  Spectral  Coefficients 


Subroutine  rad_TS _ 

! Calculate  Source  Terms 

net_emis (imech, : )  =  spect_emis_m (imech, : , i, j ) - 
spect_absb_m (imech, : , i, j ) *spect_int (imech, : ,  dir) 

call  simpson (Q_rad_dir (imech, i, j ) , wave_length_m ( : , i, j ) , net_emis (imech, : ) , nwa ve, ier) 
Q_rad_dir (imech, i, j )  =  Q_rad_dir (imech, i, j ) *1 . Oe-4 

! Calculate  Vibrational  Source  Terms 

! net_emis_vib (imech, :)  =  alpha_vib_m (imech, :, i, j ) *net_emis (imech, :) 

!  call 

simpson (Q_vib_dir (imech, i, j ) , wave_length_m ( : , i, j ) , net_emis_vib (imech, : ) , nwave, ier) 

! Q_vib_dir (imech, i, j )  =  Q_vib_dir (imech, i, j ) *1 . Oe-4 

if (dir . eq . 1 ) then 

! Calculate  Total  Emisssion  Coefficient 
call 

simpson (tot_emis (imech, i, j ) , wave_length_m ( : , i, j ) , spect_emis_m (imech, : , i, j ) , nwave, ier) 
tot_emis (imech, i, j ) =tot_emis (imech, i, j ) *le-4 

! Calculate  Total  Absorption  Coefficient 
call 

simpson (tot_absb (imech, i, j ) , wave_length_m ( : , i, j ) , spect_absb_m (imech, : , i, j ) , nwave, ier) 
tot_absb (imech, i, j ) =tot_absb (imech, i, j ) *le-4 


Table  21.  Calculate  Total  Intensities 


Subroutine  radJTS 


ICacluate  Total  Intensities 
spect_int (1,  :  ,  dir) =0 . 0 
norm_int ( 1 , i , j , dir) =0 . 0 
do  imech=2 , nmech 


norm  int  ( 1 , i 
spect  int(l, 
end  do 

j,dir)=norm  int (1 , i, j , dir) +norm  int (imech, i, j , dir) 

,dir)=spect  int (1, : , dir) +spect  int (imech, :, dir) 

Table  22.  Calculate  Source  Terms 

Subroutine  rad_TS 

do  idir  =  l,ndir 

Q  rad  tot  (:,:,: 
Q  rad  vib  (:,:,: 

=  Q  rad  tot  (:,:,: )  +  4*pi/real (ndir) *Q  rad  dir (:,:,: ) 

=  Q  rad  vib (:,:,: )  +  4*pi/real (ndir) *Q  vib  dir (:,:,: ) 

end  do 
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Table  23.  Basic  Outline  of  Finite  Volume  Method  Radiation  Solver;  important  aspects  of 
the  subroutine  rad_FVM  are  further  described  in  the  tables  indicated  below. 


Subroutine  md_FVM_ _ 

subroutine  rad_FVM (PARAMETERS , X, Y, T, Tv_l , Tv_2 , . . . , Tel, nums,molwt) 

-Declare  variables 

-Initialize  variables 

-Calculate  emis,  absb 

do  imech=l , nmech 
do  idir  =  l,ndir 

-Calculate  LHS  Matrix  and  RHS  vector  (See  Table  24) 

-Enforce  BCs  (See  Table  25) 

-Solve  Linear  System  (See  Tables  26  and  27) 

end  do  !idir 
end  do  ! imech 

-Calculate  source  terms 
end  subroutine 
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Table  24.  Calculate  =view  factor4  Matrix  (LHS)  and  Source  Term  Vector  (RHS) 


Subroutine  rad_FVM 


do  i_region=2 , 1 , -1 

do  i=2 , ILD 
do  j  =1 ,  JLD 

if (i_region . eq. 1) then 

i j  =int (  ( i-2  )  * JLD  +  j  +  (ILD-1)*JLD  ) 

else 

i j  =int (  ( ILD-i ) * JLD  +  j  ) 

end  if 

do  kdir=l , 4 

!  ***CALCULATE  LHS  MATRIX  (DIAGONAL  ELEMENTS) 

if (dot_prod (i, j , idir, kdir, i_region) . it . 0 . 0)  then 

SELECT  CASE (kdir) 

CASE (1) 

F_PENTA(1,  ij  )  =F_PENTA(1,  ij  )  -abs  (dot_prod  ( i ,  j  ,  idir ,  1 ,  i_region)  ) 

CASE (2) 

E_PENTA  (1,  ij  )  =E_PENTA  (1,  ij  )  -abs  (dot_prod  (i,  j  ,  idir,  2,  i_region)  ) 

CASE (3) 

C_PENTA(1, ij ) =C_PENTA(1, ij ) -abs (dot_prod (i, j , idir, 3, i_region) ) 

CASE  (4) 

A_PENTA  (1,  ij  )  =A_PENTA  (1,  ij  )  -abs  (dot_prod  (i,  j  ,  idir,  4,  i_region)  ) 

END  SELECT 

else 

D_PENTA  (1,  ij  )  =D_PENTA  (1,  ij  )  Tabs  (dot_prod  (i,  j  ,  idir,  kdir,  i_region)  ) 
end  if 

end  do  !kdir 
if (idir . eq. 3) then 

D_PENTA  (1,  ij  )  =D_PENTA  (1,  ij  )  +2 . 0D0*abs  (dot_prod  (i,  j  ,  idir,  5,  i_region)  ) 
end  if 


!  ***CALCULATE  RHS  VECTOR 

! if (imech.eq. 1) then 

! B_PENTA ( 1 , i j )  =  spect_emis_m (imech, : , i , j ) *vol (i , j ) *n_dir (idir, 4 ) 

!  else 

! B_PENTA ( 1 , i j )  =  spect_emis_m (imech, : , i , j ) *vol (i, j ) *n_dir (idir , 4 ) -& 

! &  spect_int (1, : , i, j , idir) *spect_absb_m (imech, : , i, j ) *vol (i, j ) *n_dir (idir, 4) 
lend  if 


if (imech.eq. 1) then 

B_PENTA ( 1 , i j )  =  tot_emis (imech, i, j ) *vol (i, j ) *n_dir (idir, 4 ) 
else 

B_PENTA ( 1 , i j )  =  tot_emis (imech, i, j ) *vol (i, j ) *n_dir (idir, 4 ) -  & 

&  norm_int (1, i, j , idir) *tot_absb (imech, i, j ) *vol (i, j ) *n_dir (idir, 4) 

end  if 


j  ***FINISH  LHS  MATRIX 

D_PENTA (1, ij ) =D_PENTA (1, ij ) +tot_absb (imech, i, j ) *vol (i, j ) *n_dir (idir, 4) 

end  do  ! j 
end  do  !  i 


end  do  ! i  region 


Table  25.  Enforce  Boundary  Conditions 


Subroutine  rad_FVM 

! Outflow 

if (i . eq. ILD  . 

and.  i  region. eq 

1 ) then 

B  PENTA ( 1 ,  i  j  ) 

=  B  PENTA (l,ij) 

D  PENTA (l,ij) 

=  D  PENTA (l,ij) 

F  PENTA (l,ij) 

=  0.0D0 

C  PENTA (l,ij) 

=  C  PENTA (l,ij) 

A  PENTA (l,ij) 

=  A  PENTA  (1,  ij  ) 

E_PENTA(1,  ij  ) 

=  E_PENTA(l,ij) 

! Outflow 

if (i . eq. ILD  . 

and.  i  region. eq 

2 ) then 

B  PENTA (l,ij) 

=  B  PENTA (l,ij) 

D  PENTA (l,ij) 

=  D  PENTA (l,ij) 

F  PENTA (l,ij) 

=  F  PENTA (l,ij) 

C  PENTA (l,ij) 

=  C  PENTA (l,ij) 

A  PENTA ( 1 ,  i  j  ) 

=  A  PENTA  (1,  ij  ) 

E_PENTA ( 1 , i j ) 

=  0.0D0 

!  Wall 

if ( j . eq. 1) then 

B  PENTA (l,ij) 

=  B  PENTA (l,ij) 

D  PENTA (l,ij) 

=  D  PENTA (l,ij) 

F  PENTA ( 1 ,  i  j  ) 

=  F  PENTA (l,ij) 

C  PENTA (l,ij) 

=  C  PENTA (l,ij) 

A  PENTA (l,ij) 

=  0.0D0 

E_PENTA(1, ij ) 
end  if 

=  E_PENTA(l,ij) 

! Inflow 

if ( j . eq. JLD) then 

B  PENTA (l,ij) 

=  B  PENTA (l,ij) 

D  PENTA (l,ij) 

=  D  PENTA (l,ij) 

F  PENTA ( 1 , i j ) 

=  F  PENTA (l,ij) 

C  PENTA (l,ij) 

=  0.0D0! 

A  PENTA (l,ij) 

=  A  PENTA (l,ij) 

E  PENTA ( 1 ,  i  j  ) 
end  if 

=  E  PENTA (l,ij) 
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Table  26.  Call  Linear  Solver;  Respecify  Solution  m  Terms  of  Global  Discretization 

Indices 


Subroutine  rad_FVM _ 

call  band  (B_PENTA,  E_PENTA,  A_PENTA,  D_PENTA,  C_PENTA, 

&  F_PENTA, 2*int ( JLD) +1 , npenta, nmech, nwave) 

do  j  =1 ,  JLD 
do  i=2 , ILD 

i j  =int (  ( i-2 ) * JLD  +  j  +  (ILD-1)*JLD  ) 

norm_int (imech, i, j , idir)  =  B_PENTA ( 1 , i j ) 

if (idir . eq. 3) then 

i j  =int (  ( i-2 ) * JLD  +  j  +  (ILD-1)*JLD  ) 

norm_int (imech, i, j , 3)  =  B_PENTA ( 1 , i j ) 

i j  =int (  ( ILD-i ) * JLD  +  j  ) 

norm_int (imech, i, j , 4 )  =  B_PENTA ( 1 , i j ) 

end  if 


norm_int (imech, 1 ,  j , idir )  =  norm_int (imech, 2 , j , idir) 
norm_int (imech, 1 ,  j ,  4 )  =  norm_int (imech, 2 ,  j ,  4 ) 


Table  27.  Linear  Solver  Used  for  Implicit  FVMR  Scheme 

Subroutine  band _ 

subroutine  band (B,E,A,D,C,F,M,N, nmech, nwave) 
implicit  none 

integer  ::  nmech, nwave, iwave,  ia,  ib 
real*8  ::  E  (N)  ,  A  (N)  ,  D  (N)  ,  C  (N)  ,  F  (N)  ,  & 

! &  B (nwave, N) 

&  B  (N) 

integer  g, h, i , j , k, m, n, r 
real*8  : :  aa (n, m) 
real*8  : :  eps 

write(*,*)  'in  band' 

nwave  =  1  ! added  for  FVM  mod 
iwave  =  1  ! added  for  FVM  mod 


r  =  (m+1 ) / 2 


eps 

=  1.0D-: 

aa=0 . 0D0 

aa  ( 

,1  ) 

=E 

aa  ( 

,r-l) 

=A 

aa  ( 

,r  ) 

=D 

aa  ( 

r  r+1) 

=C 

aa  ( 

,m  ) 

=F 

do  z 

>0  k  = 

1, 

if (  abs  (aa(k,r))  .le.  eps) 
aa(k,r)  =  1 . ODO/aa (k, r) 
h  =  r-1 
i  =  k+1 

If (h . it . 1  .or.  j  . gt .  n)  goto  20 


aa (i, h) 


=  aa (i, h) *aa (k, r) 
j  =  h+1 
g  =  r+1 

if (g.gt.m  .or.  j . gt .  (r+n-i)  )  goto  40 

aa(i,j)  =  aa(i,j)  -  aa (i, h) *aa (k, g) 

j  =  j+1 

g  =  g+1 
goto  30 

continue 
i=i  +  l 
h=h-l 
goto  10 

continue 

do  1=1 , n 

write(69,*)  (aa (i,  j  ) ,  j=r-l,r+l) 
end  do 

Forward  Elimination 


do  100  k  =  l,n-l 
i=k+l 
j=r-l 

if ( j . It . 1  .or.  i.gt.n)  goto  100 
! do  iwave  =  1 , nwave 

!b (iwave, i)  =  b (iwave, i)  -  aa (i, j ) *b (iwave, k) 
b  (i)  =  b  (i)  -  aa  (i,  j  )  *b  (k) 

! end  do 
i  =  i+1 
j  =  j-1 

goto  110 
continue 


Back  Substitution 

do  120  k  =  n, 1 , -1 
i=k+l 
j=r+l 

if(j.gt.m  .or.  i.gt.n)  goto  140 
! do  iwave  =  1, nwave 

!b (iwave, k)  =  b (iwave, k)  -  aa (k, j ) *b (iwave, i) 
! end  do 

b  (k)  =  b(k)  -  aa(k,  j)*b(i) 

1— i+1 
j=j+l 
goto  130 
continue 

! do  iwave  =  1, nwave 
!b (iwave, k)  =  b (iwave, k) *aa (k, r) 
b  (k)  =  b  (k)  *aa  (k,  r) 

! if (abs (b (iwave, k) )  .it.  1.0D-14)  b (iwave, k) =0 . 0D0 
! end  do 
continue 


end  subroutine  band 


Table  28.  Calculate  Source  Terms  for  Flow  Solver 


Subroutine  rad_couple _ 

subroutine  rad_couple ( ILD, JLD) 

use  rad_parameters 
use  rad  vars 


! Q_rad_s (1,  i,  j  )  =  total  E  rad  source  term 
! Q_rad_s (2 , i , j )  =  02  Evib  rad  source  term 

! Q_rad_s (3, i, j )  =  N2  Evib  rad  source  term 

! Q_rad_s ( 4 , i ,  j )  =  NO  Evib  rad  source  term 

!  Calculate  Source  Terms 
do  iatoms=l,18 

if (atom_rads ( 1 ,  iatoms ) . ne . '  '  ) then 

do  imech=l,nmech 
! Calculate  source  term  for  E 

if (atom_rads (1, iatoms) . eq.mech_name (imech) ) then 
Q_rad_s  (1,  : ,  : )  =  Q_rad_s  (1,  : ,  : )  +  Q_rad_tot  (imech,  :  ,  : ) 
end  if 
end  do 

end  if 
end  do 

do  iatoms=l,40 

if (diatom_bands ( 1 , iatoms ) .ne. 1  '  ) then 

do  imech=l , nmech 
! Calculate  source  term  for  E 

if (diatom_bands (1, iatoms) . eq.mech_name (imech) ) then 
Q_rad_s (1,  : ,  : )  =  Q_rad_s (1, : , : )  +  Q_rad_tot (imech, : , : ) 

! Calculate  source  term  for  Evib_02 
if ( (diatom_bands (1, iatoms) . eq. '02  ' ) .and. 

(diatom_bands (2 , iatoms) .ne. ' cont ' ) ) then 

Q_rad_s (2, : , : )  =  Q_rad_s (2, : , : )  +  Q_rad_vib (imech, : , : ) 
end  if 

! Calculate  source  term  for  Evib_N2 
if ( (diatom_bands (1, iatoms) . eq. 'N2  ' ) . and. 

(diatom_bands (2 , iatoms ) . ne . ' cont ' ) ) then 

Q_rad_s (3, : , : )  =  Q_rad_s (3, : , : )  +  Q_rad_vib (imech, : , : ) 


! Calculate  source  term  for  Evib_N0 
if ( (diatom_bands ( 1 , iatoms ) . eq . ' NO  ' ) . and . 
(diatom_bands (2 , iatoms) .ne. ' cont ' ) ) then 

Q_rad_s (4, : , : )  =  Q_rad_s (4, : , : )  +  Q_rad_vib (imech, : , : ) 
end  if 

end  if 
end  do 

end  if 
end  do 

write (*,260)  maxval (Q_rad_s ) 

260  format ( 'maxval (Q_rad_s) = ' , IpelO . 3, ' W/m3 ' ) 

end  subroutine 


Table  29.  Assign  Local  Variables  to  Shared  Variables 


Subroutine  mdipac _ 

call  emis_absb  (See  Table  30) 

spect_emis (1, : ) =spect%emis 
spect_absb (1, : ) =spect%absb 
do  isp=l , num_diatoms 

if (diatoms (isp) %name . eq. 'N2 ' )  Ni_N2  =  diatoms (isp) %state_pop 
if (diatoms (isp) %name . eq.  ' 02 '  )  Ni_02  =  diatoms (isp) %state_pop 
if (diatoms (isp) %name . eq. 'NO ' )  Ni_N0  =  diatoms (isp) %state_pop 
end  do 

do  isp=l,num_atoms 

if (atoms (isp) %name . eq. ' N ' )  Ni_N  =  atoms (isp) %state_pop 
if (atoms (isp) %name . eq. ' 0 ' )  Ni_0  =  atoms (isp) %state_pop 
end  do 


Table  30.  Calculate  Species  Contributions  to  Spectral  Coefficients 


Subroutine  emfcjjibsb;  indicates  modifications  to  baseline  code. _ 

!  bound-bound  radiation;  atomic  (See  Table  32) 
ibb  =  0 
do  i=l , 18 

if ( (atoms (isp) %name . eq . atom_rads ( 1 , i ) ) . and . (atom_rads (2 , i ) . eq . ' bb ' ) ) then 
ibb  =  1 

!  write (*,*)  ' ibb= ' , ibb 

do  j=l,nmech 

if (mech_name ( j ) . eq . atom_rads ( 1 , i ) )  then 
imech= j 

!  write ( * , * )  ' imech= '  ,  j ,  ' mech_name= ' , mech_name ( j ) 

end  if 

end  do 

end  if 
end  do 

if(ibb.eq.l)  then 

if (atoms (isp) %name . ne .' H  ')  call  atom_bb(isp,  atoms,  t,  spect , imech) 
if (atoms (isp) %name . eq. ' H  ')  call  h_bb(isp,  atoms,  t,  spect) 
end  if 

!  bound-bound  radiation;  diatomic  (See  Table  33) 

ibb  4  0 
do  1=1,40 

if ( (diatoms (isp) %name . eq . diatom_bands ( 1 , i ) ) . and . (diatom_bands (2 , i ) . ne . ' cont ' ) ) then 
ibb  =  1 

do  j=l,nmech 

if (mech_name  ( j )  . eq . diatom_bands ( 2 , i ) )  then 
imech=  j 

!  write  ( *  ,  * )  '  imech= '  ,  imech ,  '  mech_name= '  , mech_name  ( imech) 

end  if 

end  do 

end  if 
end  do 

if(ibb.eq.l)  call  diatom_bb (isp,  diatoms,  dens,  t,  spect,  diatom_bands , imech) 
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Table  3 1 .  Assign  Nonequilibrium  Temperatures  and  Number  Densities 


Subroutine  emisjibsb _ 

! Assign  Temps  an  Num  Densities;  diatomic  sp 

if (diatoms (isp) %name . eq. 'N2  ')  then 

t (isp) %vib%val  =  tvib_N2 (inode) 

diatoms (isp) %dens_diatom  =  1.0e-6  *  concN2 (inode) 

diatoms (isp) %dens_atoml  =  1.0e-6  *  concN (inode) ;  diatoms (isp) %dens_atom2  =  & 
&  1.0e-6  *  concN (inode) 

endif 

if (diatoms (isp) %name . eq. ' NO  ')  then 

diatoms (isp) %dens_diatom  =  1.0e-6  *  concNO (inode) 

t (isp) %vib%val  =  tvib_NO (inode) 

diatoms (isp) %dens_atoml  =  1.0e-6  *  concN (inode) ;  diatoms (isp) %dens_atom2  =  & 
&  1.0e-6  *  concO (inode) 

endif 

if (diatoms (isp) %name . eq. ' 02  ')  then 

t (isp) %vib%val  =  tvib_02 (inode) 

diatoms (isp) %dens_diatom  =  1.0e-6  *  conc02 (inode) 

diatoms (isp) %dens_atoml  =  1.0e-6  *  concO (inode) ;  diatoms (isp) %dens_atom2  =  & 
&  1.0e-6  *  concO (inode) 

endif 


Table  32.  Calculate  Bound-Bound  Radiation  from  Atomic  Species 


Subroutine  atomJ)b;  indicates  modifications  to  baseline  code. _ 

i *****Variable  belonging  to  common  block  /rad/****************************************** 
integer  : :  imech,  nmech 
character ( 4 )  : :  mech_name ( 68 ) 

real *8 ,  pointer  : :  spect_emis (:,:), spect_absb (:,:), alpha_vib ( : , : ) 
real*8,  pointer  ::  wave_length ( : ) 

common/ rad/ spect_emis ,  spect_absb , wave_length , nmech , me ch_name , alpha_vib 

I*************************************************************************************** 

do  m  =  nstart,nend 

dlam  =  1 . OdO/ (1 . OdO/spect%min  -  spect%itv  *  (ncentr  -  l.OdO)) 

& 

&  -  1 . OdO/ ( 1 . OdO/spect%min  -  spect%itv  *  (m  -  l.OdO)) 

emission  =  e  * ( ( ( 1 . OdO-widthl/widthv)  *  exp (-2 . 772* (dlam/widthv) **2 ) +  & 

&  (widthl/widthv) /( 1 . OdO+4 . OdO* (dlam/widthv) **2 . OdO )  +  0 . 01 6d0* (widthl/widthv) * 

& 

&  (1 . OdO-widthl/widthv) * (exp (-0 . 4d0* (abs (dlam) /widthv) **2 . 25d0)  -  10. OdO/ 

& 

&  (10 . 0d0+ (abs (dlam) /widthv) **2 . 25d0) ) ) /denom) 

spect%emis (m)  =  spect%emis (m)  +  emission 

blam  =  1.1904d-16  *  ax/ ( ( 1 . 0d-8*spect%wavel (m) ) **5 . OdO* ( 1 . OdO  -  ax)) 
spect%absb (m)  =  spect%absb (m)  +  emission/blam 

spect_emis (imech, m) =emission 
spect_absb (imech, m) =emission/blam 

end  do  ! over  m 
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Table  33.  Calculate  Bound-Bound  Radiation  from  Diatomic  Species 


Subroutine  diatomjjb;  indicates  modifications  to  baseline  code. 


I *****Variable  belonging  to  common  block  /rad/****************************************** 
integer  : :  imech,  nmech 
character (4)  ::  mech_name (68) 

real *8 ,  pointer  : :  spect_emis (:,:), spect_absb (:,:), alpha_vib ( : , : ) 
real*8,  pointer  ::  wave_length ( : ) 

common/ rad/ spect_emis  ,  spect_absb ,  wave_length ,  nmech ,  me ch_name ,  alpha_vib 

!  Sigma (upper) -Sigma (lower)  transition 

if ( (diatoms (isp) %di_lev (  diatoms (isp) %di_line (tr) %up_state  ) %lambda . eq . 0 )  & 

&  . and. (diatoms (isp) %di_lev (  diatoms (isp) %di_line (tr) %lo_state  )% lambda  & 

&  . eq. 0) )  then 

!  r-branch  ( j  +l-> j ) 

do  j  =  l,maxj 

s_j j  =  real(j,  prec) 

call  calc_diatomic_bb (isp,  diatoms, & 

&  dev,  bvu,  bvl,  dvu,  dvl,  geu,  teu,  evu,  qtot,  rel,  & 

&  j,  s_ j j ,  ’ r’,  spect,  wavelx,  emisj,  ncentr) 

j 

!  duplicate  line  profile  of  band  origin 

do  m  =  -nspred, nspred 

if ( (ncentr ( j ) +m. gt . 0) . and. (ncentr ( j ) +m. le . spect%wave_num) )  then 
emission  =  emisj  (j)  *  y (m) 

ax  =  exp(-  1 . 43877*1 . 0e8/ (spect%wavel (ncentr (j ) +m) *  & 

&  t (isp) %el_tr%val) ) 

blam  =  1.1904e-16  *  ax/ ( (1 . 0e-8*spect%wavel (ncentr (j ) +m) ) **5*  & 

&  ( 1 . 0  -  ax) ) 

absorption  =  emission/blam 

alpha_vib (imech, m)  =  exp (-1 . 43877*evu/t (isp) %vib%val) /qtot 

spect_emis (imech, ncentr ( j) +m) =emission 
spect_absb ( imech , ncentr ( j ) +m) =absorption 

spect%emis (ncentr (j )  +  m)  =  spect%emis (ncentr (j )  +  m)  +  emission 
spect%absb (ncentr (j )  +  m)  =  spect%absb (ncentr (j )  +  m)  +  absorption 
end  if 
enddo 
end  do 

This  modification  repeated  for  other  transition  types  (i.e.,  S-P,  P-S,  P-P). _ 
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Appendix  B:  State  Populations 


Distance  from  Wall  (cm)  Distance  from  Wall  (cm) 


Figure  20.  a)  State  Populations  of  N,  1634.0  seconds  (left  uncoupled,  right  coupled) 


Distance  from  Wall  (cm)  Distance  from  Wall  (cm) 


Figure  20.  b)  State  Populations  of  O,  1634.0  seconds  (left  uncoupled,  right  coupled) 


Distance  from  Wall  (cm)  Distance  from  Wall  (cm) 


Figure  20.  c)  State  Populations  of  N2,  1634.0  seconds  (left  uncoupled,  right  coupled) 


140 


Figure  20.  d)  State  Populations  of  O2,  1634.0  seconds  (left  uncoupled,  right  coupled) 


Figure  20.  e)  State  Populations  of  NO,  1634.0  seconds  (left  uncoupled,  right  coupled) 


Distance  from  Wall  (cm)  Distance  from  Wall  (cm) 


Figure  21.  a)  State  Populations  of  N,  1636.0  seconds  (left  uncoupled,  right  coupled) 
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Distance  from  Wall  (cm)  Distance  from  Wall  (cm) 

Figure  21.  b)  State  Populations  of  O,  1636.0  seconds  (left  uncoupled,  right  coupled) 


Distance  from  Wall  (cm)  Distance  from  Wall  (cm) 


Figure  21.  c)  State  Populations  of  N2,  1636.0  seconds  (left  uncoupled,  right  coupled) 


Distance  from  Wall  (cm) 


Distance  from  Wall  (cm) 


Figure  21.  d)  State  Populations  of  O2,  1636.0  seconds  (left  uncoupled,  right  coupled) 
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Distance  from  Wall  (cm) 


Figure  21.  e)  State  Populations  of  NO,  1636.0  seconds  (left  uncoupled,  right  coupled) 


Distance  from  Wall  (cm)  Distance  from  Wall  (cm) 

Figure  22.  a)  State  Populations  of  N,  1640.5  seconds  ( left  uncoupled,  right  coupled) 


Distance  from  Wall  (cm)  Distance  from  Wall  (cm) 


Figure  22.  b)  State  Populations  of  O,  1640.5  seconds  (left  uncoupled,  right  coupled) 
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— I — X1SIGg+ 
— 9 —  A3SIGU+ 
— # —  B3Plg 
x  W3DELU 
— h —  Bp3SIGu- 
— 0 —  apISIGu- 
— A — alPIg 


wIDELu 

Ap5SIGg+ 

G3DELg 


Figure  22.  c)  State  Populations  of  N2,  1640.5  seconds  (left  uncoupled,  right  coupled) 
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Figure  22. 


d)  State  Populations  of  02,  1640.5  s  (left  uncoupled,  right  coupled) 
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Figure  22. 


e)  State  Populations  of  NO,  1640.5  s  (left  uncoupled,  right  coupled) 
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